Week 2: Introduction to Likelihood Inference
PS 818 - Data Analysis with Statistical Models
September 14, 2026
Likelihood Inference
\[ \DeclareMathOperator*{\argmin}{arg\,min} \DeclareMathOperator*{\argmax}{arg\,max} \]
Previously
- Linear regression
- “Agnostic” regression and the Best Linear Predictor
- “Classical” OLS regression theory
- Key intuitions
- Regression as a projection of \(Y\) into the column space of \(\mathbf{X}\)
- Frisch-Waugh-Lovell theorem
- Throughout, we assumed as little as possible about the distribution of \(Y\).
This week
- What is a statistical model?
- What is a likelihood?
- Properties of the maximum likelihood estimator
- How do we actually compute the MLE – and its standard errors?
- This last part turns out to be much less automatic than it looks.
Statistical models
Defining a statistical model
- In the regression setting we tried to make as few assumptions about the data-generating process as possible.
- Our goal is just to estimate and conduct inference on \(E[Y|X]\).
- But what if we wanted to make further probabilistic statements about other quantities beyond \(\beta\)?
- (e.g.) Can we provide a distribution for \(Y_{n+1}\), the “next” observation given \(X_{n+1}\)?
- If we’re willing to make more assumptions about the data-generating process, we can do a lot more!
- Statistical models specify the data-generating process in terms of systematic and stochastic components.
- Systematic elements are functions of known constants and unknown parameters
- Stochastic elements are draws from probability distributions
- We will be primarily working with parametric models
- The data are assumed to come from a particular family of probability distributions
- The “structure” of the model is fixed (the number of parameters does not grow with the size of the data).
The linear model
It is common to see the linear model written in its fully parametric form.
Stochastic:
\[y_i \sim \text{Normal}(\mu_i, \sigma^2)\]
Systematic:
\[\mu_i = x_i^{\prime}\beta\]
What’s assumed to be known?
- \(\mathbf{x}\)
What’s assumed to have a particular distribution?
- \(y\)
We are interested in estimating and conducting inference on the parameters: \(\beta\) and (less importantly) \(\sigma^2\).
General model notation
We can specify a broad set of models for \(y_i\) using this framework
Stochastic
\[y_i \sim f(\theta_i, \alpha)\]
Systematic
\[\theta_i = g(x_i, \beta)\]
General model notation
What are these quantities?
- \(y_i\) is a random variable
- \(f()\) denotes the distribution of that random variable
- \(\theta_i\) and \(\alpha\) are parameters of that distribution
- \(g()\) is some function
- \(x_i\) are observed, known constants (e.g. regressors)
- \(\beta\) are parameters of interest
Later, we will spend some time with a particular class of models called “Generalized Linear Models” where the systematic component has the form
\[\theta_i = g(x_i^{\prime}\beta)\]
Types of distributions
- Normal
Continuous on an unbounded support \((-\infty, \infty)\)
Two parameters: Mean \(\mu\) and Variance \(\sigma^2\)
\[f(x;\mu, \sigma) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left\{- \frac{1}{2}\left(\frac{x - \mu}{\sigma}\right)^2 \right\}\]
Types of distributions
- Poisson
Discrete, defined on the support of the natural numbers (positive integers + zero)
Single parameter: Mean and Variance \(\lambda\)
\[f(x; \lambda) = \frac{\lambda^x\exp\{-\lambda\}}{x!}\]
Types of distributions
- Binomial
- Discrete, defined on the support of integers from \(\{0, 1, 2, \dotsc, n\}\) (models the sum of repeated i.i.d. coin flips)
- Two parameters: \(p\) probability of success in \(n\) trials
- Special case where \(n=1\) trial is typically called the “Bernoulli”
Likelihood inference
Learning about the unknown
- Suppose that I want to learn about an unobserved parameter from a sample of observations.
- I assume a particular statistical model for the data
- Example: The linear model from earlier!
- Frequentist approach
- Unobserved parameters are fixed constants
- Data are random variables
- We want to construct an estimator that is a function of the data and which has desirable properties
- Consistency: As our sample size gets larger, the estimator converges (in probability) to the target parameter.
- Asymptotic normality: In large samples, the distribution of our estimator is normal (ideally with a variance we can estimate as well!)
- Can we come up with a generic framework that yields a “good” estimator for a large class of statistical models?
The Likelihood Function
Consider data \(\mathbf{y} = \{y_1, y_2, \dotsc, y_n\}\) which we assume comes from a known distribution with density \(f(\mathbf{y} | \theta)\) with unknown parameter \(\theta\)
- \(f()\) could be Bernoulli, Normal, Poisson, Gamma, Beta, etc… – in this setting it is a known distribution
- \(\theta\) is an unknown parameter that lies in some space \(\Theta\) of possible values.
The likelihood function is a function of \(\theta\) that is evaluated at the observed values of \(\mathbf{y}\)
\[\mathcal{L}(\theta| \mathbf{y}) = f(\mathbf{y} | \theta)\]
Given some input \(\theta\), the likelihood function returns the density of the observed data evaluated at that value.
The likelihood function is not a probability density
- a pdf is a function that takes \(x\) as an input
- a likelihood is a function that takes \(\theta\) as an input
The Likelihood Function
The likelihood function is not \(f(\theta | \mathbf{y})\)
- That statement doesn’t even make sense in a frequentist framework – parameters are constants
Even when we (later) move to a Bayesian framework, \(f(\theta | \mathbf{y}) \neq f(\mathbf{y} | \theta)\)
Remember Bayes’ Rule:
\[f(\theta | \mathbf{y}) = \frac{f(\mathbf{y} | \theta)}{f(\mathbf{y})} \times f(\theta)\]
Or alternatively, since \(f(\mathbf{y})\) is just a normalizing constant, you’ll see this written as:
\[\underbrace{f(\theta | \mathbf{y})}_{\text{posterior}} \propto \underbrace{f(\mathbf{y}|\theta)}_{\text{likelihood}} \times \underbrace{f(\theta)}_{\text{prior}}\]
The Likelihood Function
We often make the assumption that our data are independently and identically distributed
- \(y_i \underset{{\text{i.i.d}}}{\sim} f(y_i | \theta)\)
This allows us to factor the likelihood
\[\mathcal{L}(\theta|\mathbf{y}) = f(\mathbf{y} | \theta) = \prod_{i=1}^n f(y_i | \theta)\]
Since our eventual goal will be to find an optimum of the likelihood, we can apply any strictly increasing function to it, since that preserves the location of the maximum.
- The main one we’ll apply is the logarithm
- Why? Because logs turn annoying-to-work-with products into easier-to-work-with sums!
Notation: sums vs. averages
For theory reasons, you’ll see two forms of the log-likelihood
The total log-likelihood
\[\ell_n(\theta) = \sum_{i=1}^n \log f(y_i | \theta)\]
The average log-likelihood
\[\bar{\ell}_n(\theta) = \frac{1}{n}\ell_n(\theta) = \frac{1}{n}\sum_{i=1}^n \log f(y_i | \theta)\]
These have the same maximizer
- But the average lets us make law-of-large-numbers arguments - this will matter later on!
Notation: Score and Hessian
The score for a single observation is the gradient of its log-density with respect to \(\theta\)
\[s_i(\theta) = \frac{\partial}{\partial\theta} \log f(y_i | \theta)\]
For a \(k\)-dimensional \(\theta\), \(s_i(\theta)\) is a \(k \times 1\) vector.
The Hessian for a single observation is the matrix of second derivatives
\[H_i(\theta) = \frac{\partial^2}{\partial\theta\partial\theta^{\prime}} \log f(y_i | \theta)\]
Notation: Score and Hessian
Summing and averaging as before:
\[S_n(\theta) = \sum_{i=1}^n s_i(\theta) \qquad \bar{S}_n(\theta) = \frac{1}{n}\sum_{i=1}^n s_i(\theta)\]
\[H_n(\theta) = \sum_{i=1}^n H_i(\theta) \qquad \bar{H}_n(\theta) = \frac{1}{n}\sum_{i=1}^n H_i(\theta)\]
\(S_n(\theta)\) is exactly the gradient of \(\ell_n(\theta)\), and \(H_n(\theta)\) its Hessian.
MLE
We want to come up with an estimator \(\hat{\theta}\) for the parameter \(\theta\)
- \(\hat{\theta}\) is a function of the data (like a sample mean or OLS coefficient)
- Is there a principled way to pick \(\hat{\theta}\) that has provably “good” properties across a wide variety of models?
The Maximum Likelihood Estimator (MLE) \(\hat{\theta}\) is defined as the value of \(\theta\) that maximizes the log-likelihood
\[\hat{\theta} = \argmax_{\theta \in \Theta} \ \ell_n(\theta)\]
Under i.i.d. observations (and regularity conditions), the MLE has some desirable properties:
- Consistent: \(\hat{\theta} \overset{p}{\to} \theta_0\), the true parameter
- Asymptotically normal: \(\sqrt{n}(\hat{\theta} - \theta_0) \overset{d}{\to} \text{Normal}\big(0, \mathcal{I}(\theta_0)^{-1}\big)\)
- Asymptotically efficient: its variance attains the Cramér-Rao lower bound
The score has mean zero
First, we’ll start with some properties of the score.
One property that comes up frequently: evaluated at the parameter that generated the data, the score has mean zero
\[\mathbb{E}_{\theta}\big[s_i(\theta)\big] = 0\]
Start from the definition of the score
\[\mathbb{E}_{\theta}\big[s_i(\theta)\big] = \int \left[\frac{\partial}{\partial\theta} \log f(y|\theta)\right] f(y|\theta) \, dy\]
Apply the chain rule to \(\log f\), then cancel \(f(y|\theta)\)
\[= \int \frac{\partial}{\partial\theta} f(y|\theta) \, dy\]
The score has mean zero
Swap the derivative and the integral (regularity condition - we assume we’re able to do this)
\[= \frac{\partial}{\partial\theta} \int f(y|\theta) \, dy\]
The density integrates to one
\[= \frac{\partial}{\partial\theta} \, 1 \; = \; 0\]
Consistency
By the weak law of large numbers, for each \(\theta\)
\[\bar{\ell}_n(\theta) = \frac{1}{n} \sum_{i=1}^n \log f(y_i | \theta) \overset{p}{\rightarrow} \mathbb{E}_{\theta_0}[\log f(y_i | \theta)] \]
The MLE maximizes the left-hand side, so (under regularity conditions) \(\hat{\theta}\) converges to whatever maximizes \(\mathbb{E}_{\theta_0}[\log f(y_i | \theta)]\)
We can show that the true parameter \(\theta_0\) maximizes \(\mathbb{E}_{\theta_0}[\log f(y_i | \theta)]\)
Consistency
Consider any other value of \(\theta\) and apply properties of logs
\[\mathbb{E}_{\theta_0}[\log f(y_i | \theta)] - \mathbb{E}_{\theta_0}[\log f(y_i | \theta_0)] = \mathbb{E}_{\theta_0}\left[\log \frac{f(y_i|\theta)}{f(y_i|\theta_0)}\right]\]
By Jensen’s inequality (\(\log\) is concave),
\[\mathbb{E}_{\theta_0}\left[\log \frac{f(y_i|\theta)}{f(y_i|\theta_0)}\right] \le \log \mathbb{E}_{\theta_0}\left[\frac{f(y_i|\theta)}{f(y_i|\theta_0)}\right] = \log \int f(y|\theta) \, dy = \log 1 = 0\]
So \(\mathbb{E}_{\theta_0}[\log f(y_i | \theta)] \le \mathbb{E}_{\theta_0}[\log f(y_i | \theta_0)]\) for all \(\theta\)
- And we get equality at \(\theta = \theta_0\): \(\theta_0\) maximizes \(\mathbb{E}_{\theta_0}[\log f(y_i | \theta)]\)
Therefore the MLE is consistent: \(\hat{\theta} \overset{p}{\to} \theta_0\).
Consistency
Note a few important conditions for identification
Identifiability
- No “plateaus” in the log-likelihood.
- Two different values of \(\theta\) can’t both maximize the log-likelihood
\[f(\cdot\,|\theta) \neq f(\cdot\,| \theta_0) \ \forall \ \theta \neq \theta_0\]
- This is exactly what makes the inequality above strict aside from \(\theta = \theta_0\)
Fixed parameter space
- The dimensionality of \(\theta\) stays fixed and does not depend on \(n\)
- Violated (e.g.) in ideal point models (new legislators mean new ideal points)
\(\theta_0\) is in the interior of the parameter space
- Otherwise the score need not be zero at the optimum.
Information
After consistency, we want the variance and asymptotic distribution of \(\hat{\theta}\).
- To do this, we need a quantity called the Fisher information.
Define the per-observation information as the variance of the score at the true parameter \(\theta_0\)
\[\mathcal{I}(\theta_0) = \mathbb{E}_{\theta_0}\big[s_i(\theta_0) s_i(\theta_0)^{\prime}\big]\]
Under some further regularity conditions this equals the negative expected Hessian
\[\mathcal{I}(\theta_0) = -\mathbb{E}_{\theta_0}\big[H_i(\theta_0)\big]\]
Intuitively: captures the curvature of the log-likelihood around the truth.
- Sharper peak \(\Rightarrow\) more information \(\Rightarrow\) smaller variance.
Information
The total (sample) information for \(n\) i.i.d. observations is
\[\mathcal{I}_n(\theta_0) = -\mathbb{E}_{\theta_0}\big[H_n(\theta_0)\big] = n \, \mathcal{I}(\theta_0)\]
Asymptotic distribution
Now we can show that the MLE is asymptotically normal.
\[\sqrt{n}(\hat{\theta} - \theta_0) \overset{d}{\rightarrow} \text{Normal}\big(0, \mathcal{I}(\theta_0)^{-1}\big)\]
Our first step is to take a first-order Taylor approximation of the score around the true parameter \(\theta_0\), written as a function of our estimator \(\hat{\theta}\)
\[\bar{S}_n(\hat{\theta}) \approx \bar{S}_n(\theta_0) + \bar{H}_n(\theta_0)(\hat{\theta} - \theta_0)\]
Using the fact that the score is zero at the MLE
\[0 = \bar{S}_n(\theta_0) + \bar{H}_n(\theta_0)(\hat{\theta} - \theta_0)\]
Rearranging gives
\[(\hat{\theta} - \theta_0) = \big[-\bar{H}_n(\theta_0)\big]^{-1}\bar{S}_n(\theta_0)\]
Asymptotic distribution
Multiplying by \(\sqrt{n}\)
\[\sqrt{n}(\hat{\theta} - \theta_0) = \big[-\bar{H}_n(\theta_0)\big]^{-1}\big[\sqrt{n}\,\bar{S}_n(\theta_0)\big]\]
First, by the law of large numbers applied to the average Hessian,
\[-\bar{H}_n(\theta_0) = -\frac{1}{n} \sum_{i=1}^n H_i(\theta_0) \overset{p}{\to} -\mathbb{E}_{\theta_0}\big[H_i(\theta_0)\big] = \mathcal{I}(\theta_0)\]
Second, by the central limit theorem applied to the mean-zero scores,
\[\sqrt{n} \, \bar{S}_n(\theta_0) = \frac{1}{\sqrt{n}} \sum_{i=1}^n s_i(\theta_0) \overset{d}{\to} \text{Normal}\big(0, \, \text{Var}_{\theta_0}(s_i(\theta_0))\big) = \text{Normal}\big(0, \mathcal{I}(\theta_0)\big)\]
Asymptotic distribution
Lastly, by Slutsky’s theorem,
\[\sqrt{n}(\hat{\theta} - \theta_0) \overset{d}{\to} \text{Normal}\big(0, \mathcal{I}(\theta_0)^{-1}\mathcal{I}(\theta_0)\mathcal{I}(\theta_0)^{-1}\big) = \text{Normal}\big(0, \mathcal{I}(\theta_0)^{-1}\big)\]
So our large-sample approximation for the distribution of the MLE is
\[\hat{\theta} \overset{a}{\sim} \text{Normal}\big(\theta_0, \ \mathcal{I}_n(\theta_0)^{-1}\big)\]
Standard errors are the square roots of the diagonal of \(\mathcal{I}_n(\hat{\theta})^{-1}\).
Can another estimator do better asymptotically?
- It turns out the answer is no
- The MLE asymptotically attains the Cramér-Rao lower bound.
Cramér-Rao lower bound (CRLB)
An important result in statistics relates to the smallest possible variance that an unbiased estimator could reach
- That is, if we’re not willing to accept more bias, how much precision can we get out of the data?
- This ends up being related to the information!
The Cramér-Rao lower bound on the variance of any unbiased estimator \(\hat{\theta}\) of \(\theta\) is
\[\text{Var}(\hat{\theta}) \succeq \mathcal{I}_n(\theta)^{-1}\]
We’ll focus the proof on the scalar case
- For vector \(\theta\), “\(\succeq\)” means the difference \(\text{Var}(\hat{\theta}) - \mathcal{I}_n(\theta)^{-1}\) is positive semi-definite.
Cramér-Rao lower bound (CRLB)
Under unbiasedness
\[\mathbb{E}_{\theta}[\hat{\theta}(\mathbf{y}) - \theta] = \int(\hat{\theta}(\mathbf{y}) - \theta) f(\mathbf{y}| \theta)d\mathbf{y} = 0\]
Since unbiasedness holds for every value of \(\theta\), differentiating with respect to \(\theta\) also gives zero
\[\frac{\partial}{\partial\theta} \int(\hat{\theta}(\mathbf{y}) - \theta) f(\mathbf{y} |\theta)d\mathbf{y} = 0\]
Under the regularity conditions that let us exchange differentiation and integration
\[\int \frac{\partial}{\partial\theta} \bigg[(\hat{\theta}(\mathbf{y}) - \theta) f(\mathbf{y}| \theta)\bigg]d\mathbf{y} = 0\]
Cramér-Rao lower bound (CRLB)
Product rule
\[\int (\hat{\theta}(\mathbf{y}) - \theta) \frac{\partial}{\partial\theta} f(\mathbf{y}| \theta) \, d\mathbf{y} - \int f(\mathbf{y}| \theta)d\mathbf{y} = 0\]
PDFs integrate to 1 over their support
\[\int (\hat{\theta}(\mathbf{y}) - \theta) \frac{\partial}{\partial\theta} f(\mathbf{y}| \theta) \, d\mathbf{y} = 1\]
The same chain rule identity as before
\[\frac{\partial}{\partial \theta} f(\mathbf{y}| \theta) = f(\mathbf{y}| \theta) \frac{\partial}{\partial \theta} \log f(\mathbf{y}|\theta) = f(\mathbf{y}|\theta) \, S_n(\theta)\]
Cramér-Rao lower bound (CRLB)
Plugging back in,
\[\int (\hat{\theta}(\mathbf{y}) - \theta) \, S_n(\theta) \, f(\mathbf{y}| \theta)d\mathbf{y} = 1\]
By the “law of the unconscious statistician” this is just an expectation
\[\mathbb{E}_{\theta}\big[(\hat{\theta}(\mathbf{y}) - \theta) \, S_n(\theta)\big] = 1\]
The Cauchy-Schwarz inequality for expectations says that for random variables \(U\) and \(V\)
\[\big|\mathbb{E}[UV]\big|^2 \le \mathbb{E}[U^2]\,\mathbb{E}[V^2]\]
Apply it with \(U = \hat{\theta}(\mathbf{y}) - \theta\) and \(V = S_n(\theta)\)
\[1 = \big|\mathbb{E}_\theta[U V]\big|^2 \le \mathbb{E}_{\theta}\big[(\hat{\theta}(\mathbf{y}) - \theta)^2\big] \times \mathbb{E}_{\theta}\big[S_n(\theta)^2\big]\]
Cramér-Rao lower bound (CRLB)
What are the pieces here?
- \(\mathbb{E}_{\theta}[(\hat{\theta}(\mathbf{y}) - \theta)^2]\) is the mean squared error. Because we assumed \(\hat{\theta}\) is unbiased, MSE is the variance.
- \(\mathbb{E}_{\theta}[S_n(\theta)^2]\) is the total information \(\mathcal{I}_n(\theta)\)
So
\[\text{Var}(\hat{\theta}) \times \mathcal{I}_n(\theta) \ge 1 \quad \Longrightarrow \quad \text{Var}(\hat{\theta}) \ge \frac{1}{\mathcal{I}_n(\theta)}\]
Compare this to the large-sample distribution of the MLE from earlier:
\[\hat{\theta}_{MLE} \overset{a}{\sim} \text{Normal}\big(\theta_0, \ \mathcal{I}_n(\theta_0)^{-1}\big)\]
The MLE attains the bound asymptotically.
The Normal model
Deriving the log-likelihood
Let’s consider the classical normal regression model with homoskedastic errors
\[y_i \sim \text{Normal}(\mu_i, \sigma^2), \qquad \mu_i = x_i^{\prime}\beta\]
Under i.i.d. the joint density factors, and logs turn the product into a sum
\[\ell_n(\beta, \sigma^2) = \log \prod_{i=1}^n f(y_i| \beta, \sigma^2) = \sum_{i=1}^n \log f(y_i| \beta, \sigma^2)\]
Plug in the Normal density for each observation
\[\ell_n(\beta, \sigma^2) = \sum_{i=1}^n \log \bigg[\frac{1}{\sigma\sqrt{2\pi}} \exp\left\{- \frac{1}{2}\left(\frac{y_i - x_i^\prime\beta}{\sigma}\right)^2 \right\}\bigg]\]
Properties of logs – note \(\log(\sigma\sqrt{2\pi}) = \log\sigma + \tfrac{1}{2}\log(2\pi)\)
\[\ell_n(\beta, \sigma^2) = \sum_{i=1}^n \left[-\log(\sigma) - \tfrac{1}{2}\log(2\pi) -\frac{1}{2}\bigg(\frac{y_i - x_i^\prime\beta}{\sigma}\bigg)^2\right]\]
Deriving the log-likelihood
The \(\tfrac{1}{2}\log(2\pi)\) term does not depend on \(\beta\) or \(\sigma\), so it cannot change the location of the maximum. Ignoring it,
\[\ell_n(\beta, \sigma^2) = -\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - x_i^\prime\beta)^2 -n\log(\sigma)\]
And the corresponding average log-likelihood
\[\bar{\ell}_n(\beta, \sigma^2) = \frac{1}{n} \sum_{i=1}^n \bigg[-\frac{1}{2\sigma^2} (y_i - x_i^\prime\beta)^2 - \log(\sigma)\bigg]\]
What’s familiar here?
- For fixed \(\sigma\), maximizing \(\ell_n\) over \(\beta\) is the same as minimizing the sum of squared residuals.
- So the MLE for \(\beta\) is the OLS estimator!
Deriving the score
Let’s derive the score. The per-observation log-density, again ignoring the constant, with residuals \(e_i = y_i - x_i^\prime\beta\):
\[\log f(y_i|\beta, \sigma) = -\frac{e_i^2}{2\sigma^2} - \log(\sigma)\]
Differentiating with respect to \(\beta\) (chain rule, \(\partial e_i / \partial\beta = -x_i\)):
\[\frac{\partial}{\partial\beta} \log f(y_i|\beta,\sigma) = \frac{x_i \, e_i}{\sigma^2}\]
Differentiating with respect to \(\sigma\):
\[\frac{\partial}{\partial\sigma} \log f(y_i|\beta,\sigma) = \frac{e_i^2}{\sigma^3} - \frac{1}{\sigma}\]
Setting the sum of the first to zero gives the OLS normal equations \(\mathbf{X}^\prime\mathbf{e} = 0\)
- The second gives \(\hat{\sigma}^2 = \frac{1}{n}\sum_i e_i^2\). Note the MLE is not the unbiased estimator of the \(\sigma^2\)
Deriving the Hessian
Differentiating once more gives the blocks of \(H_n\):
\[\frac{\partial^2 \ell_n}{\partial\beta\partial\beta^\prime} = -\frac{\mathbf{X}^\prime\mathbf{X}}{\sigma^2}, \qquad \frac{\partial^2 \ell_n}{\partial\beta \, \partial\sigma} = -\frac{2\mathbf{X}^\prime\mathbf{e}}{\sigma^3}, \qquad \frac{\partial^2 \ell_n}{\partial\sigma^2} = \frac{n}{\sigma^2} - \frac{3\sum_i e_i^2}{\sigma^4}\]
Evaluated at the MLE
- \(\mathbf{X}^\prime\mathbf{e} = 0\), so the cross-derivative is exactly zero – \(\hat\beta\) and \(\hat\sigma\) are orthogonal.
- \(\sum_i e_i^2 = n\hat\sigma^2\), so the \(\sigma\sigma\) entry collapses to \(-2n/\hat\sigma^2\).
So the exact Hessian at the MLE is block diagonal:
\[H_n(\hat{\theta}) = \begin{bmatrix} -\mathbf{X}^\prime\mathbf{X}/\hat\sigma^2 & 0 \\ 0 & -2n/\hat\sigma^2 \end{bmatrix}\]
Which gives us closed-form standard errors
\[\widehat{\text{Var}}(\hat\beta) = \hat\sigma^2(\mathbf{X}^\prime\mathbf{X})^{-1}\]
Our familiar regression variance-covariance matrix for \(\hat{\beta}\)!
Fitting the model
- How does Clinton’s 2016 county-level vote share vary by the share of residents who are non-Hispanic white?
Fitting the model
- We’ll use a numerical optimizer implemented in R
- The
maxLikpackage is a nice wrapper for the various optimizers built in to R
- The
- The main routine,
maxLik, takes a function whose first argument is the parameter vector.- We’ll initially give it only the log-likelihood, so every derivative it needs is computed numerically.
- First, write a function returning the log-likelihood evaluated at each observation
- For optimization we reparameterize \(\sigma = \exp(\eta)\), so the optimizer can range over all of \(\mathbb{R}\) and never has to worry about \(\sigma \le 0\).
Fitting the model
- We initialize our \(\mathbf{X}\) matrix, our outcome \(Y\) and some starting values for the optimizer.
- Then pass the function
norm_loglikand all other arguments tomaxLik().- We’ll choose
method = "NR"to select the Newton-Raphson optimizer
- We’ll choose
[1] 67.799 -0.462 2.505- Compare the point estimates to OLS:
Standard errors
- We use a plug-in estimator for the standard errors.
- Calculate the Hessian at the MLE and plug it into our asymptotic variance formula we derived earlier (invert the negative hessian).
- Compare these to what we get from conventional OLS standard errors
term se_maxLik se_ols pct_error
(Intercept) (Intercept) 0.49869 0.9059 -44.9
pctWhiteNH pctWhiteNH 0.00657 0.0112 -41.3- These should be identical (aside from degrees of freedom correction in
lm) - Instead the slope standard error is off by about 41%.
The answer isn’t stable!
- Try some different starting values (vary \(\log\sigma\)’s initial point)
se_from_start <- function(ls0, grad=NULL){
fit <- maxLik(norm_loglik, grad=grad, start=c(0, 0, ls0), method="NR",
Y=Y_vec, X=X_matrix)
se <- suppressWarnings(sqrt(diag(solve(-maxLik::hessian(fit)))))
c(beta1_hat = coef(fit)[2], SE_beta1 = se[2])
}
starts <- c(2.0, 2.303, exp(1), 2.708, 3.0)
round(t(sapply(starts, se_from_start)), 6) beta1_hat SE_beta1
[1,] -0.462 0.00844
[2,] -0.462 0.00845
[3,] -0.462 0.00845
[4,] -0.462 NaN
[5,] -0.462 0.00657- The standard error swings wildly, and in some runs the negative Hessian is not even positive definite
- The truth (from OLS) is 0.011.
- The optimizer is behaving poorly
Numerical vs. analytical derivatives
- The Newton-Raphson algorithm uses the gradient and the Hessian to find the optimum iteratively.
- It can get them two ways:
- Finite differences (“numerical differentiation”)
- \(\frac{\partial f}{\partial\theta} \approx \frac{f(\theta + h) - f(\theta - h)}{2h}\) (done one parameter at a time)
- Needs only the ability to evaluate \(f\), so it works with just the likelihood as input.
- Approximate, the error depends on software-selected \(h\).
- Exact derivatives
- Just do the calculus and pass the optimizer a function that returns the score.
- We’ve already done this for the linear model!
Including the analytic gradient
- We derived the score earlier. Write this as a function that returns an \(N \times (k+1)\) matrix – one gradient row per observation.
Analytic gradient: the standard errors
term numeric analytic ols_truth
(Intercept) (Intercept) 0.49869 0.9056 0.9059
pctWhiteNH pctWhiteNH 0.00657 0.0112 0.0112- Supplying the gradient recovers the OLS standard errors essentially exactly.
maxLikis still finite-differencing to get the Hessian – but now it differences the exact gradient - less room for numerical error.
And now it is stable
- Across different starting values
beta1_hat SE_beta1
[1,] -0.462 0.0112
[2,] -0.462 0.0112
[3,] -0.462 0.0112
[4,] -0.462 0.0112
[5,] -0.462 0.0112- The same answer from every start.
Standardizing \(Y\) and \(X\)
- Another solution is to make the inputs better behaved
- Normalize \(Y\) and \(X\) by subtracting the mean and dividing by the standard deviation.
- We’ll compare the two approaches
- The discrepancy disappears
se_z_numeric <- sqrt(diag(solve(-maxLik::hessian(mle_z_numeric))))
se_z_analytic <- sqrt(diag(solve(-maxLik::hessian(mle_z_analytic))))
data.frame(
term = c("(Intercept)", "pctWhiteNH", "log(sigma)"),
numeric = se_z_numeric,
analytic = se_z_analytic,
pct_diff = 100*(se_z_numeric - se_z_analytic)/se_z_analytic
) term numeric analytic pct_diff
1 (Intercept) 0.0144 0.0144 -0.000879
2 pctWhiteNH 0.0144 0.0144 -0.001603
3 log(sigma) 0.0127 0.0127 0.001887Getting back to interpretable units
To get back to the units of interest, undo the re-scaling:
\[\hat\beta_1^{raw} = \hat\beta_1^{z} \times \frac{s_Y}{s_X}, \qquad \widehat{\text{SE}}(\hat\beta_1^{raw}) = \widehat{\text{SE}}(\hat\beta_1^{z}) \times \frac{s_Y}{s_X}\]
quantity from_z ols_truth
Estimate slope -0.4617 -0.4617
Std. Error SE(slope) 0.0112 0.0112“Robust” standard errors
Relaxing the information matrix equality
Return to the asymptotic variance derivation. We had the sandwich
\[\text{Var}\big(\sqrt{n}(\hat{\theta}-\theta_0)\big) \approx \underbrace{\big[-\mathbb{E}_{\theta_0}[H_i(\theta_0)]\big]^{-1}}_{\text{bread}} \underbrace{\text{Var}_{\theta_0}\big(s_i(\theta_0)\big)}_{\text{meat}} \underbrace{\big[-\mathbb{E}_{\theta_0}[H_i(\theta_0)]\big]^{-1}}_{\text{bread}}\]
Previously we used the information matrix equality – \(\text{meat} = \text{bread}^{-1} = \mathcal{I}(\theta_0)\) – to collapse this to \(\mathcal{I}(\theta_0)^{-1}\).
- That equality holds only if the model is correctly specified. If it isn’t, the two differ, and only the sandwich is consistent.
Instead, we can estimate each piece separately with sample analogues:
\[\widehat{\text{Var}}(\hat{\theta}) = \bigg[-\sum_{i=1}^n H_i(\hat\theta)\bigg]^{-1}\bigg[\sum_{i=1}^n s_i(\hat{\theta}) s_i(\hat{\theta})^{\prime}\bigg]\bigg[-\sum_{i=1}^n H_i(\hat\theta)\bigg]^{-1}\]
Consistent under misspecification of the outcome distribution, but finite-sample properties can be poor absent corrections (e.g. \(\frac{n}{n-k}\) scaling known as HC1).
Implementing the sandwich
sandwich::estfunreturns the per-observation gradient
[1] 3114 3[1] -1.0e-06 -5.1e-05 3.1e-05- Assembling the three pieces:
term classical robust robust_hc1
1 (Intercept) 0.9056 1.0348 1.0351
2 pctWhiteNH 0.0112 0.0122 0.0122
3 log(sigma) 0.0127 0.0138 0.0139Conclusion
- Statistical models
- Describe the data-generating process in terms of systematic and stochastic components
- In a parametric model, we assume the data \(\mathbf{y}\) come from a known distribution with unknown parameters \(\theta\)
- Likelihood
- A function of \(\theta\) evaluated at the observed data \(\mathbf{y}\), equal to \(f(\mathbf{y}|\theta)\)
- How relatively plausible are my observed results in a world where the true DGP parameter were \(\theta\)?
- Not \(f(\theta|\mathbf{y})\). That quantity only makes sense in the Bayesian context and requires a prior \(f(\theta)\)
- Maximum Likelihood Estimator
- A technique for constructing a good estimator whenever we can write down a well-behaved likelihood
- Consistency, asymptotic normality, asymptotic efficiency – all large sample properties
Next week
- More on maximum likelihood estimators
- Invariance, the delta method, likelihood ratio / Wald / score tests
- How do optimizers work
- Newton’s method, gradient ascent/descent
- Generalized linear models
- Link functions and the exponential family
- Binary outcome models
- Logit and probit as likelihood models
- Readings: Chapters 13-14 of Gelman, Hill, and Vehtari (2020)
References
PS 818 - University of Wisconsin-Madison