n_obs <- 13
xbar_obs <- 72 / 13
se_rate <- sqrt(xbar_obs / n_obs)
z_wald <- (xbar_obs - 4) / se_rate
c(
lower = xbar_obs - qnorm(0.975) * se_rate,
upper = xbar_obs + qnorm(0.975) * se_rate,
z = z_wald,
p_value = 2 * pnorm(-abs(z_wald))
)
#> lower upper z p_value
#> 4.2591657 6.8177574 2.3570226 0.0184221Introduction to Maximum Likelihood Inference
1 Overview of maximum likelihood estimation
These notes are derived primarily from (Dobson and Barnett 2018, chaps. 1–5), with some material from (McLachlan and Krishnan 2007) and (Casella and Berger 2002).
1.1 The likelihood function
NoteNotation for the likelihood function
Sources write the likelihood function in several ways, all meaning the same function:
- \(\mathcal{L}(\theta)\);
- \(\mathcal{L}(\tilde{x}; \theta)\) or \(\mathcal{L}(\theta; \tilde{x})\);
- \(\mathcal{L}_{\tilde{x}}(\theta)\) or \(\mathcal{L}_{\theta}(\tilde{x})\);
- \(\mathcal{L}(\tilde{x}\mid \theta)\).
These notes mostly write \(\mathcal{L}(\theta)\), leaving the data implicit, to emphasize that the likelihood is a function of the parameters, with the data held fixed at their observed values.
ProofProof
Proof. \[ \begin{aligned} \mathcal{L}(\theta) &\stackrel{\text{def}}{=}\operatorname{p}(X_1 = x_1, \ldots, X_n = x_n \mid \theta) && \text{(definition of likelihood)}\\ &= \prod_{i=1}^n \operatorname{p}(X_i = x_i \mid \theta) && \text{(definition of mutual independence)} \end{aligned} \]
ProofProof
Proof. By Theorem 1, \(\mathcal{L}(\theta) = \prod_{i=1}^n \operatorname{p}(X_i = x_i \mid \theta)\), and by Definition 3, each factor \(\operatorname{p}(X_i = x_i \mid \theta)\) is \(\mathcal{L}_i(\theta)\).
1.2 The maximum likelihood estimate
1.3 Finding the maximum of a function
From calculus: if \(f(x)\) is differentiable, its maximum over an interval of input values can occur only at an endpoint of the interval or at a critical point. At a critical point \(x_0\), \(f''(x_0) < 0\) is sufficient for \(x_0\) to be a local maximum, but not necessary: \(f(x) = -x^4\) has a maximum at \(x_0 = 0\), where \(f''(0) = 0\). For a function of a vector, a negative definite Hessian matrix at a critical point is sufficient for a local maximum.
1.4 Directly maximizing the likelihood function for independent data
To find the maximizer of the likelihood function, we solve \(\mathcal{L}'(\theta) = 0\) for \(\theta\). For mutually independent data, Equation 1 gives:
\[ \begin{aligned} \mathcal{L}'(\theta) &= \frac{\partial}{\partial \theta} \mathcal{L}(\theta)\\ &= \frac{\partial}{\partial \theta} \prod_{i=1}^n \operatorname{p}(X_i = x_i \mid \theta) \end{aligned} \tag{3}\]
Equation 3 is the derivative of a product of \(n\) factors, which takes \(n - 1\) applications of the product rule and produces \(n\) terms. The log-likelihood avoids this work.
1.5 The log-likelihood function
ProofProof
Proof. The natural logarithm is strictly increasing on \((0, \infty)\), so for any \(\theta_1\) and \(\theta_2\), \(\mathcal{L}(\theta_1) \ge \mathcal{L}(\theta_2)\) if and only if \(\log \mathcal{L}(\theta_1) \ge \log \mathcal{L}(\theta_2)\), that is, \(\ell(\theta_1) \ge \ell(\theta_2)\). So \(\theta^*\) satisfies \(\mathcal{L}(\theta^*) \ge \mathcal{L}(\theta)\) for every \(\theta\) if and only if it satisfies \(\ell(\theta^*) \ge \ell(\theta)\) for every \(\theta\).
ProofProof
Proof. \[ \begin{aligned} \ell(\theta) &\stackrel{\text{def}}{=}\log{\mathcal{L}(\theta)} && \text{(definition of log-likelihood)}\\ &= \log{\prod_{i=1}^n \operatorname{p}(X_i = x_i \mid \theta)} && \text{(likelihood of an independent sample)}\\ &= \sum_{i=1}^n \log{\operatorname{p}(X_i = x_i \mid \theta)} && \text{(log of a product is a sum of logs)} \end{aligned} \]
With a common distribution, \(\operatorname{p}(X_i = x_i \mid \theta) = \operatorname{p}(X = x_i \mid \theta)\).
ProofProof
Proof. \[ \begin{aligned} \ell'(\theta) &= \frac{\partial}{\partial \theta} \ell(\theta) && \text{(notation for the derivative)}\\ &= \frac{\partial}{\partial \theta} \sum_{i=1}^n \log{\operatorname{p}(X = x_i \mid \theta)} && \text{(log-likelihood of an $\operatorname{iid}$ sample)}\\ &= \sum_{i=1}^n \frac{\partial}{\partial \theta} \log{\operatorname{p}(X = x_i \mid \theta)} && \text{(derivative of a sum is the sum of derivatives)} \end{aligned} \]
Unlike Equation 3, each term involves only one observation, so no product rule is needed.
1.6 The score function
We often omit the arguments \(\tilde{x}\) and \(\theta\), writing \(\ell' \stackrel{\text{def}}{=}\ell'(\tilde{x}\mid \theta) \stackrel{\text{def}}{=}\ell'(\theta)\). Some sources write \(U\) or \(S\) for the score function instead of \(\ell'\); for example, Dobson and Barnett (2018) write \(U\). These notes use \(\ell'\), which keeps \(U\) and \(S\) free for other uses and needs no extra symbol to memorize.
In all four examples (Exercise 3, Exercise 4, Exercise 5, and Exercise 6), the score function with respect to the mean turned out to be:
\[\ell'= \frac{x - \operatorname{E}\mathopen{}\left[X\right]\mathclose{}}{\operatorname{Var}\mathopen{}\left(X\right)\mathclose{}}\]
This pattern is no coincidence. With the mean as the parameter, each of these four models is a one-parameter natural (or linear) exponential family, whose log-density is linear in \(x\); for every such family, the score with respect to the mean is \((x - \operatorname{E}\mathopen{}\left[X\right]\mathclose{})/\operatorname{Var}\mathopen{}\left(X\right)\mathclose{}\). Other members of the broader exponential family, such as the Weibull distribution with known shape \(k \ne 1\), do not have this form. Exponential-family distributions share many special properties (Hogg et al. 2019, sec. 6.7; Dobson and Barnett 2018, chap. 3).
1.7 Information matrices
The Hessian is named after the mathematician Otto Hesse.
ProofProof
Proof. \(\frac{\partial}{\partial \tilde{\theta}^{\top}}\ell\) is the \(1 \times p\) row vector whose \(j\)th entry is \(\frac{\partial}{\partial \theta_j}\ell\). Differentiating each entry with respect to the \(p \times 1\) vector \(\tilde{\theta}\) gives a \(p \times 1\) column of derivatives per entry, so \(\frac{\partial}{\partial \tilde{\theta}}\frac{\partial}{\partial \tilde{\theta}^{\top}}\ell\) is \(p \times p\), with \(ij\)th entry \(\frac{\partial}{\partial \theta_i}\frac{\partial}{\partial \theta_j}\ell\).
ProofProof
Proof. By Definition 7, \(\ell'= \frac{\partial}{\partial \tilde{\theta}}\ell\), so \(\mathopen{}\left(\ell'\right)\mathclose{}^{\top} = \frac{\partial}{\partial \tilde{\theta}^{\top}}\ell\), and \(\frac{\partial}{\partial \tilde{\theta}}\mathopen{}\left(\ell'\right)\mathclose{}^{\top} = \frac{\partial}{\partial \tilde{\theta}}\frac{\partial}{\partial \tilde{\theta}^{\top}}\ell= \ell''\) by Definition 9.
ProofProof
Proof. We give the proof for a scalar \(\theta\) and a continuous \(\tilde{X}\) with density \(\operatorname{p}(\tilde{x}\mid \theta)\); the vector case applies the same steps to each pair of entries, and a discrete \(\tilde{X}\) replaces integrals with sums.
For Equation 12:
\[ \begin{aligned} \operatorname{E}\mathopen{}\left[\ell'\right]\mathclose{} &= \int \mathopen{}\left(\frac{\partial}{\partial \theta} \log \operatorname{p}(\tilde{x}\mid \theta)\right)\mathclose{} \operatorname{p}(\tilde{x}\mid \theta) \, d\tilde{x} && \text{(definition of expectation)}\\ &= \int \frac{\frac{\partial}{\partial \theta} \operatorname{p}(\tilde{x}\mid \theta)}{\operatorname{p}(\tilde{x}\mid \theta)} \operatorname{p}(\tilde{x}\mid \theta) \, d\tilde{x} && \text{(chain rule for $\log$)}\\ &= \int \frac{\partial}{\partial \theta} \operatorname{p}(\tilde{x}\mid \theta) \, d\tilde{x} && \text{(cancel $\operatorname{p}(\tilde{x}\mid \theta)$)}\\ &= \frac{\partial}{\partial \theta} \int \operatorname{p}(\tilde{x}\mid \theta) \, d\tilde{x} && \text{(exchange derivative and integral)}\\ &= \frac{\partial}{\partial \theta} 1 && \text{(a density integrates to 1)}\\ &= 0 \end{aligned} \]
For Equation 13, differentiate the identity \(0 = \int \mathopen{}\left(\frac{\partial}{\partial \theta} \log \operatorname{p}(\tilde{x}\mid \theta)\right)\mathclose{} \operatorname{p}(\tilde{x}\mid \theta) \, d\tilde{x}\) from the first line with respect to \(\theta\):
\[ \begin{aligned} 0 &= \int \frac{\partial}{\partial \theta}\mathopen{}\left[\mathopen{}\left(\frac{\partial}{\partial \theta} \log \operatorname{p}(\tilde{x}\mid \theta)\right)\mathclose{} \operatorname{p}(\tilde{x}\mid \theta)\right]\mathclose{} d\tilde{x} && \text{(exchange derivative and integral)}\\ &= \int \mathopen{}\left(\frac{\partial^2}{\partial \theta^2} \log \operatorname{p}(\tilde{x}\mid \theta)\right)\mathclose{} \operatorname{p}(\tilde{x}\mid \theta) \, d\tilde{x} + \int \mathopen{}\left(\frac{\partial}{\partial \theta} \log \operatorname{p}(\tilde{x}\mid \theta)\right)\mathclose{} \frac{\partial}{\partial \theta}\operatorname{p}(\tilde{x}\mid \theta) \, d\tilde{x} && \text{(product rule)}\\ &= \operatorname{E}\mathopen{}\left[\ell''\right]\mathclose{} + \int \mathopen{}\left(\frac{\partial}{\partial \theta} \log \operatorname{p}(\tilde{x}\mid \theta)\right)\mathclose{}^2 \operatorname{p}(\tilde{x}\mid \theta) \, d\tilde{x} && \text{($\frac{\partial}{\partial \theta}\operatorname{p}= \mathopen{}\left(\frac{\partial}{\partial \theta}\log \operatorname{p}\right)\mathclose{} \operatorname{p}$)}\\ &= \operatorname{E}\mathopen{}\left[\ell''\right]\mathclose{} + \operatorname{E}\mathopen{}\left[\ell'^2\right]\mathclose{} && \text{(definition of expectation)} \end{aligned} \]
So \(\mathcal{I}(\theta) = \operatorname{E}\mathopen{}\left[-\ell''\right]\mathclose{} = \operatorname{E}\mathopen{}\left[\ell'^2\right]\mathclose{}\), and because \(\operatorname{E}\mathopen{}\left[\ell'\right]\mathclose{} = 0\), \(\operatorname{E}\mathopen{}\left[\ell'^2\right]\mathclose{} = \operatorname{E}\mathopen{}\left[\ell'^2\right]\mathclose{} - \mathopen{}\left(\operatorname{E}\mathopen{}\left[\ell'\right]\mathclose{}\right)^2\mathclose{} = \operatorname{Var}\mathopen{}\left(\ell'\right)\mathclose{}\).
Sources disagree on the symbols for the observed and expected information (Table 1).
1.8 Asymptotic distribution of the maximum likelihood estimate
ProofProof
Proof. The proof is beyond the scope of these notes; see (Lehmann 1999, Theorem 7.5.2, p. 501) and (Newey and McFadden 1994).
These conditions guarantee a consistent root of the score equation; that root is the global maximizer of the likelihood under further conditions, for example when the log-likelihood is strictly concave, as for Poisson data, whose Hessian is negative for every \({\lambda}\) (Example 5).
Theorem 9 involves the unknown \(\tilde{\theta}\), so to use it we estimate \(\mathcal{I}(\tilde{\theta})\), by either the expected information at the MLE, \(\mathcal{I}(\hat\theta_{\text{ML}})\), or the observed information at the MLE, \(I(\tilde{x}; \hat\theta_{\text{ML}})\). Either way, the estimated standard error of the \(k\)th entry of \(\hat\theta_{\text{ML}}\) is:
\[ \mathop{\widehat{\operatorname{SE}}}\nolimits\mathopen{}\left(\hat\theta_k\right)\mathclose{} = \sqrt{\mathopen{}\left[\mathopen{}\left(\hat{\mathcal{I}}\right)^{-1}\mathclose{}\right]\mathclose{}_{kk}} \]
where \(\hat{\mathcal{I}}\) is whichever estimate of \(\mathcal{I}(\tilde{\theta})\) we chose.
Using the observed information is often more convenient, and there are settings where it is provably better by some criteria (Efron and Hinkley 1978).
1.9 Quantifying uncertainty about MLEs
1.9.1 Confidence intervals for MLEs
By Theorem 9, \((\hat\theta_k - \theta_k)/\mathop{\widehat{\operatorname{SE}}}\nolimits\mathopen{}\left(\hat\theta_k\right)\mathclose{}\) has approximately a standard Gaussian distribution in large samples, so the Wald interval is an approximate confidence interval for \(\theta_k\). For a 95% interval, \(z_{0.975} \approx 1.96\).
1.9.2 Wald tests
1.9.3 Likelihood ratio tests for MLEs
ProofProof
Proof. We sketch the argument for a scalar parameter and the point null hypothesis \(H_0: \theta= \theta_0\), which imposes \(q = 1\) constraint. Then the restricted MLE is \(\hat\theta_0 = \theta_0\), and \(\Lambda = 2\mathopen{}\left(\ell(\hat\theta_{\text{ML}}) - \ell(\theta_0)\right)\mathclose{}\). For the full proof, see (Wilks 1938) or (Dobson and Barnett 2018, sec. 5.7).
Step 1: expand the log-likelihood around the MLE. The MLE maximizes \(\ell\) at an interior point, so \(\ell'(\hat\theta_{\text{ML}}) = 0\). A second-order Taylor expansion of \(\ell(\theta_0)\) around \(\hat\theta_{\text{ML}}\) gives:
\[ \begin{aligned} \ell(\theta_0) &\approx \ell(\hat\theta_{\text{ML}}) + \ell'(\hat\theta_{\text{ML}})(\theta_0 - \hat\theta_{\text{ML}}) + \frac{1}{2}\ell''(\hat\theta_{\text{ML}})(\theta_0 - \hat\theta_{\text{ML}})^2 && \text{(Taylor expansion)}\\ &= \ell(\hat\theta_{\text{ML}}) + \frac{1}{2}\ell''(\hat\theta_{\text{ML}})(\theta_0 - \hat\theta_{\text{ML}})^2 && \text{($\ell'(\hat\theta_{\text{ML}}) = 0$)} \end{aligned} \]
Rearranging:
\[ \begin{aligned} \Lambda &= 2\mathopen{}\left(\ell(\hat\theta_{\text{ML}}) - \ell(\theta_0)\right)\mathclose{} && \text{(definition of $\Lambda$)}\\ &\approx -\ell''(\hat\theta_{\text{ML}})(\hat\theta_{\text{ML}}- \theta_0)^2 && \text{(substitute the expansion)}\\ &= I(\hat\theta_{\text{ML}})(\hat\theta_{\text{ML}}- \theta_0)^2 && \text{(observed information is the negative Hessian)} \end{aligned} \]
Step 2: replace the observed information by the expected information. Under the regularity conditions, \(I(\hat\theta_{\text{ML}})/\mathcal{I}(\theta_0) \to 1\) in probability, by the law of large numbers and the consistency of \(\hat\theta_{\text{ML}}\), so:
\[\Lambda \approx \mathcal{I}(\theta_0)(\hat\theta_{\text{ML}}- \theta_0)^2\]
Step 3: recognize a squared standard Gaussian variable. Define \(Z \stackrel{\text{def}}{=}\sqrt{\mathcal{I}(\theta_0)}\,(\hat\theta_{\text{ML}}- \theta_0)\), so that \(\Lambda \approx Z^2\). By Theorem 9, \(\hat\theta_{\text{ML}}\ \dot{\sim} \ \operatorname{N}\mathopen{}\left(\theta_0, \mathcal{I}(\theta_0)^{-1}\right)\mathclose{}\) under \(H_0\), so \(Z \ \dot{\sim} \ \operatorname{N}\mathopen{}\left(0, 1\right)\mathclose{}\), and therefore \(\Lambda \approx Z^2\) has approximately a \(\chi^2_1\) distribution.
With \(q\) constraints, the same argument in matrix form makes \(\Lambda\) approximately a sum of \(q\) squared, independent standard Gaussian variables, which has a \(\chi^2_q\) distribution.
Equivalently, in terms of nested models: if a full model \(M_1\) has \(p\) free parameters, and a nested model \(M_0 \subset M_1\), obtained by imposing \(q\) constraints on \(M_1\), has \(p_0 = p - q\) free parameters, then when \(M_0\) is true, \(\Lambda = 2\mathopen{}\left(\ell_{M_1}(\hat\theta_{\text{ML}}) - \ell_{M_0}(\hat\theta_0)\right)\mathclose{}\) converges in distribution to \(\chi^2_q\).
1.9.4 Exact and approximate tests
| Inference goal | Exact test (Gaussian outcomes) | Null distribution | Approximate test (MLE) | Approximate null distribution |
|---|---|---|---|---|
| Single coefficient | t-test | \(T_{n-p}\) | Wald z-test | \(Z \sim N(0,1)\) |
| Linear combination of coefficients | t-test | \(T_{n-p}\) | Wald z-test | \(Z \sim N(0,1)\) |
| Nested models (\(q\) constraints) | Partial F-test | \(F_{q,\, n-p}\) | Likelihood ratio test or Wald test | \(\chi^2_q\) |
| One-sample mean | One-sample t-test | \(T_{n-1}\) | z-test | \(Z \sim N(0,1)\) |
| Two-sample means | Two-sample t-test (pooled variance) | \(T_{n_1+n_2-2}\) | z-test | \(Z \sim N(0,1)\) |
| \(K\)-group means (ANOVA) | F-test | \(F_{K-1,\, n-K}\) | Likelihood ratio test or Wald test | \(\chi^2_{K-1}\) |
The \(t\) and \(F\) distributions are defined in Statistical Inference. The exact tests assume \(Y_i \ \sim_{\perp\!\!\!\perp}\ N(\mu_i, \sigma^2)\), with a common variance. The approximate tests hold asymptotically, for any model that is correctly specified and satisfies the regularity conditions of Theorem 9, Gaussian or not.
1.9.5 Prediction intervals
Suppose \(X_1, \ldots, X_n \ \sim_{\operatorname{iid}}\ \operatorname{N}\mathopen{}\left(\mu, \sigma^2\right)\mathclose{}\) with \(\sigma^2\) known, and we want to predict the mean \(\bar X^*\) of \(m\) new observations from the same distribution, independent of the first \(n\). The MLE of \(\mu\) is \(\hat\mu = \bar X\), and the prediction error \(\bar X^* - \hat\mu\) is a difference of independent Gaussian variables, so:
\[ \begin{aligned} \operatorname{Var}\mathopen{}\left(\bar X^* - \hat\mu\right)\mathclose{} &= \operatorname{Var}\mathopen{}\left(\bar X^*\right)\mathclose{} + \operatorname{Var}\mathopen{}\left(\hat\mu\right)\mathclose{} && \text{(variance of a difference of independent variables)}\\ &= \frac{\sigma^2}{m} + \frac{\sigma^2}{n} && \text{(variance of a sample mean)} \end{aligned} \]
and \(\bar X^* - \hat\mu \sim \operatorname{N}\mathopen{}\left(0, \sigma^2\mathopen{}\left(\frac{1}{m} + \frac{1}{n}\right)\mathclose{}\right)\mathclose{}\). So a \(100(1 - \alpha)\%\) prediction interval for \(\bar X^*\) is:
\[\hat\mu \pm z_{1 - \alpha/2} \, \sigma \sqrt{\frac{1}{m} + \frac{1}{n}}\]
Usually \(m = 1\). The term \(1/n\) accounts for the uncertainty in \(\hat\mu\), and becomes negligible when \(n\) is much larger than \(m\).
2 Example: maximum likelihood for tropical cyclones in Australia
Adapted from (Dobson and Barnett 2018, sec. 1.6.5).
2.1 Data
Table 3 records the number of tropical cyclones in northeastern Australia during 13 November-to-April cyclone seasons, from 1956/57 to 1968/69 (Dobson and Barnett 2018, sec. 1.6.5). Figure 1 graphs the number of cyclones by season. Let \(X_i\) represent the number of cyclones in season \(i\), and \(x_i\) its observed value.
Show R code
| years | season | number |
|---|---|---|
| 1956/7 | 1 | 6 |
| 1957/8 | 2 | 5 |
| 1958/9 | 3 | 4 |
| 1959/60 | 4 | 6 |
| 1960/1 | 5 | 6 |
| 1961/2 | 6 | 3 |
| 1962/3 | 7 | 12 |
| 1963/4 | 8 | 7 |
| 1964/5 | 9 | 4 |
| 1965/6 | 10 | 2 |
| 1966/7 | 11 | 6 |
| 1967/8 | 12 | 7 |
| 1968/9 | 13 | 4 |
2.2 Exploratory analysis
Suppose we want to learn how many cyclones to expect per season.
Show R code
cyclones |>
dplyr::mutate(years = factor(years, levels = years)) |>
ggplot2::ggplot() +
ggplot2::aes(x = years, y = number, group = 1) +
ggplot2::geom_point() +
ggplot2::geom_line() +
ggplot2::xlab("Season") +
ggplot2::ylab("Number of cyclones") +
ggplot2::expand_limits(y = 0) +
ggplot2::theme(axis.text.x = ggplot2::element_text(vjust = .5, angle = 45))
Figure 1 shows no obvious trend, and no obvious correlation between adjacent seasons, so we assume that the seasons’ counts are mutually independent. We also assume that they are identically distributed, with a common distribution \(\Pr(X = x)\); the expression has no index \(i\), because the distribution is the same for every season.
Figure 2 shows the empirical distribution of the counts.
Show R code
Table 4 provides summary statistics.
Show R code
cyclones data
| seasons | total | mean | variance |
|---|---|---|---|
| 13 | 72 | 5.54 | 6.1 |
2.3 Model
We want to estimate \(\Pr(X = x)\); that is, \(\Pr(X = x)\) is our estimand.
We could estimate \(\Pr(X = x)\) for each value of \(x\) in \(0, 1, 2, \ldots\) separately (“nonparametrically”), using the fraction of our data with \(X_i = x\); but then we would be estimating an infinite set of parameters from 13 observations, and each estimate would be imprecise. A parametric model, with a few parameters, will probably do better.
2.4 Estimating the model parameters using maximum likelihood
We can estimate the parameter \(\lambda\) using maximum likelihood estimation.
SolutionSolution
Solution.
Show R code
# `cyclones` is defined earlier on the page:
# nolint next: object_usage_linter.
lik <- function(lambda, y = cyclones$number, n = length(y)) {
lambda^sum(y) * exp(-n * lambda) / prod(factorial(y))
}
lik_plot <-
ggplot2::ggplot() +
ggplot2::geom_function(fun = lik, n = 1001) +
ggplot2::xlim(min(cyclones$number), max(cyclones$number)) +
ggplot2::ylab("likelihood") +
ggplot2::xlab("lambda")
print(lik_plot)
SolutionSolution
Solution.
Show R code
# `cyclones` is defined earlier on the page:
# nolint next: object_usage_linter.
loglik <- function(lambda, y = cyclones$number, n = length(y)) {
sum(y) * log(lambda) - n * lambda - sum(lfactorial(y))
}
ll_plot <- ggplot2::ggplot() +
ggplot2::geom_function(fun = loglik, n = 1001) +
ggplot2::xlim(min(cyclones$number), max(cyclones$number)) +
ggplot2::ylab("log-likelihood") +
ggplot2::xlab("lambda")
ll_plot
2.4.1 The score function
SolutionSolution
Solution.
Show R code
# `cyclones` is defined earlier on the page:
# nolint next: object_usage_linter.
score <- function(lambda, y = cyclones$number, n = length(y)) {
(sum(y) / lambda) - n
}
ggplot2::ggplot() +
ggplot2::geom_function(fun = score, n = 1001) +
ggplot2::xlim(min(cyclones$number), max(cyclones$number)) +
ggplot2::ylab("l'(lambda)") +
ggplot2::xlab("lambda") +
ggplot2::geom_hline(yintercept = 0, col = "red")
2.4.2 The Hessian
SolutionSolution
Solution.
Show R code
# `cyclones` is defined earlier on the page:
# nolint next: object_usage_linter.
hessian <- function(lambda, y = cyclones$number, n = length(y)) {
-sum(y) / (lambda^2)
}
ggplot2::ggplot() +
ggplot2::geom_function(fun = hessian, n = 1001) +
ggplot2::xlim(min(cyclones$number), max(cyclones$number)) +
ggplot2::ylab("l''(lambda)") +
ggplot2::xlab("lambda") +
ggplot2::geom_hline(yintercept = 0, col = "red")
2.5 Finding the MLE analytically
In this example, we can find the MLE of \(\lambda\) by solving the score equation algebraically.
Call this solution of the score equation \(\tilde \lambda\) for now:
\[\tilde \lambda \stackrel{\text{def}}{=}\bar x\]
SolutionSolution
Solution.
Show R code
mle_data <- tibble::tibble(x = mle, y = loglik(mle))
ll_plot +
ggplot2::geom_point(data = mle_data, ggplot2::aes(x = x, y = y), col = "red")
Figure 8 graphs the observed information, \(I(\lambda; \tilde{x}) = -\ell''(\lambda; \tilde{x})\).
Show R code
obs_inf <- function(...) -hessian(...) # nolint: object_usage_linter.
ggplot2::ggplot() +
ggplot2::geom_function(fun = obs_inf, n = 1001) +
ggplot2::xlim(min(cyclones$number), max(cyclones$number)) +
ggplot2::ylab("I(lambda)") +
ggplot2::xlab("lambda") +
ggplot2::geom_hline(yintercept = 0, col = "red")
2.6 Finding the MLE using the Newton-Raphson algorithm
2.6.1 Iterative maximization
When we cannot solve the score equation \(\ell'(\theta) = 0\) algebraically, we can search for its solution numerically (Dobson and Barnett 2018, chap. 4).
The second line of the update uses \(I= -\ell''\) (observed information).
The update comes from approximating the score function near \({\widehat{\theta}}^*\) by its first-order Taylor polynomial:
\[ \begin{aligned} \ell'(\theta) &\approx \ell'^*(\theta)\\ &\stackrel{\text{def}}{=}\ell'({\widehat{\theta}}^*) + \ell''({\widehat{\theta}}^*)(\theta- {\widehat{\theta}}^*) \end{aligned} \]
The approximate score function \(\ell'^*(\theta)\) is linear in \(\theta\), so the approximate score equation \(\ell'^*(\theta) = 0\) is easy to solve:
\[ \begin{aligned} 0 &= \ell'({\widehat{\theta}}^*) + \ell''({\widehat{\theta}}^*)(\theta- {\widehat{\theta}}^*) && \text{(set $\ell'^*(\theta) = 0$)}\\ -\ell'({\widehat{\theta}}^*) &= \ell''({\widehat{\theta}}^*)(\theta- {\widehat{\theta}}^*) && \text{(subtract $\ell'({\widehat{\theta}}^*)$)}\\ -\mathopen{}\left(\ell''({\widehat{\theta}}^*)\right)^{-1}\mathclose{}\ell'({\widehat{\theta}}^*) &= \theta- {\widehat{\theta}}^* && \text{(multiply by $\mathopen{}\left(\ell''({\widehat{\theta}}^*)\right)^{-1}\mathclose{}$ on the left)}\\ \theta&= {\widehat{\theta}}^* - \mathopen{}\left(\ell''({\widehat{\theta}}^*)\right)^{-1}\mathclose{}\ell'({\widehat{\theta}}^*) && \text{(add ${\widehat{\theta}}^*$)} \end{aligned} \]
and the solution becomes the next guess.
The expected information is sometimes simpler to compute than the observed information.
For \(\operatorname{iid}\) data, \(\frac{1}{n}I_e(\theta; \tilde{x})\) is the sample covariance matrix (with divisor \(n\)) of the observations’ scores, so it estimates the covariance matrix of one observation’s score, \(\mathcal{I}(\theta)/n\) (Theorem 8). The empirical information needs only first derivatives, so it can be easier to compute than the observed information, and it can replace the observed information in the Newton-Raphson update.
2.6.2 Applying Newton-Raphson to the cyclone data
From Exercise 16 and Exercise 18, the score function and Hessian are:
\[ \begin{aligned} \ell'(\lambda; \tilde{x}) &= \frac{72}{\lambda} - 13\\ \ell''(\lambda; \tilde{x}) &= -\frac{72}{\lambda^2} \end{aligned} \]
So the first-order Taylor approximation of the score function around \({\widehat{\lambda}}^*\) is:
\[ \begin{aligned} \ell'(\lambda) &\approx \ell'^*(\lambda)\\ &\stackrel{\text{def}}{=}\ell'({\widehat{\lambda}}^*) + \ell''({\widehat{\lambda}}^*)(\lambda - {\widehat{\lambda}}^*)\\ &= \mathopen{}\left(\frac{72}{{\widehat{\lambda}}^*} - 13\right)\mathclose{} + \mathopen{}\left(-\frac{72}{\mathopen{}\left({\widehat{\lambda}}^*\right)^2\mathclose{}}\right)\mathclose{} (\lambda - {\widehat{\lambda}}^*) \end{aligned} \]
Figure 9 compares the score function and the approximate score function at \({\widehat{\lambda}}^*= 3\).
Show R code
# score(), hessian() and loglik() are defined earlier on the page
# nolint start: object_usage_linter.
approx_score <- function(lambda, lhat, ...) {
score(lambda = lhat, ...) +
hessian(lambda = lhat, ...) * (lambda - lhat)
}
# nolint end
point_size <- 5
plot1 <- ggplot2::ggplot() +
ggplot2::geom_function(
fun = score,
ggplot2::aes(col = "score function"),
n = 1001
) +
ggplot2::geom_function(
fun = approx_score,
ggplot2::aes(col = "approximate score function"),
n = 1001,
args = list(lhat = cur_lambda_est)
) +
ggplot2::geom_point(
size = point_size,
ggplot2::aes(
x = cur_lambda_est, y = score(lambda = cur_lambda_est),
col = "current estimate"
)
) +
ggplot2::geom_point(
size = point_size,
ggplot2::aes(x = xbar, y = 0, col = "MLE")
) +
ggplot2::xlim(min(cyclones$number), max(cyclones$number)) +
ggplot2::ylab("l'(lambda)") +
ggplot2::xlab("lambda") +
ggplot2::geom_hline(yintercept = 0)
print(plot1)
Approximating the score function by a linear function is equivalent to approximating the log-likelihood by a second-order Taylor polynomial (Figure 10):
\[ \ell^*(\lambda) \stackrel{\text{def}}{=} \ell({\widehat{\lambda}}^*) + (\lambda - {\widehat{\lambda}}^*) \ell'({\widehat{\lambda}}^*) + \frac{1}{2}\ell''({\widehat{\lambda}}^*)(\lambda - {\widehat{\lambda}}^*)^2 \]
Show R code
# nolint start: object_usage_linter.
approx_loglik <- function(lambda, lhat, ...) {
loglik(lambda = lhat, ...) +
score(lambda = lhat, ...) * (lambda - lhat) +
1 / 2 * hessian(lambda = lhat, ...) * (lambda - lhat)^2
}
# nolint end
plot_loglik <- ggplot2::ggplot() +
ggplot2::geom_function(
fun = loglik,
ggplot2::aes(col = "log-likelihood"),
n = 1001
) +
ggplot2::geom_function(
fun = approx_loglik,
ggplot2::aes(col = "approximate log-likelihood"),
n = 1001,
args = list(lhat = cur_lambda_est)
) +
ggplot2::geom_point(
size = point_size,
ggplot2::aes(
x = cur_lambda_est, y = loglik(lambda = cur_lambda_est),
col = "current estimate"
)
) +
ggplot2::geom_point(
size = point_size,
ggplot2::aes(x = xbar, y = loglik(xbar), col = "MLE")
) +
ggplot2::xlim(min(cyclones$number) - 1, max(cyclones$number)) +
ggplot2::ylab("log-likelihood") +
ggplot2::xlab("lambda")
print(plot_loglik)
Solving the approximate score equation \(\ell'^*(\lambda) = 0\) gives the next estimate:
\[ \begin{aligned} \lambda &= {\widehat{\lambda}}^*- \ell'({\widehat{\lambda}}^*) \cdot\mathopen{}\left(\ell''({\widehat{\lambda}}^*)\right)^{-1}\mathclose{}\\ &= 4.375 \end{aligned} \]
new_lambda_est <-
cur_lambda_est - score(cur_lambda_est) / hessian(cur_lambda_est)Show R code
plot2 <- plot1 +
ggplot2::geom_point(
size = point_size,
ggplot2::aes(x = new_lambda_est, y = 0, col = "new estimate")
) +
ggplot2::geom_segment(
arrow = grid::arrow(),
linewidth = 2,
alpha = .7,
ggplot2::aes(
x = cur_lambda_est,
y = approx_score(lhat = cur_lambda_est, lambda = cur_lambda_est),
xend = new_lambda_est,
yend = 0,
col = "update"
)
)
print(plot2)
We update \({\widehat{\lambda}}^*\leftarrow 4.375\) and repeat the process (Figure 12).
Show R code
plot2 +
ggplot2::geom_function(
fun = approx_score,
ggplot2::aes(col = "new approximate score function"),
n = 1001,
args = list(lhat = new_lambda_est)
) +
ggplot2::geom_point(
size = point_size,
ggplot2::aes(
x = new_lambda_est, y = score(lambda = new_lambda_est),
col = "new estimate"
)
)
We repeat this process until the log-likelihood stops changing (Table 5).
cur_lambda_est <- 3 # restart from the initial guess
tolerance <- 10^-4
max_iter <- 100
nr_info <- tibble::tibble(
iteration = 0,
lambda = cur_lambda_est,
`log(likelihood)` = loglik(cur_lambda_est),
score = score(cur_lambda_est),
hessian = hessian(cur_lambda_est)
)
for (cur_iter in 1:max_iter) {
new_lambda_est <-
cur_lambda_est - score(cur_lambda_est) / hessian(cur_lambda_est)
diff_loglik <- loglik(new_lambda_est) - loglik(cur_lambda_est)
nr_info <- nr_info |>
dplyr::bind_rows(
tibble::tibble(
iteration = cur_iter,
lambda = new_lambda_est,
`log(likelihood)` = loglik(new_lambda_est),
score = score(new_lambda_est),
hessian = hessian(new_lambda_est),
`diff(loglik)` = diff_loglik
)
)
cur_lambda_est <- new_lambda_est
if (abs(diff_loglik) < tolerance) {
break
}
}
nr_info |> knitr::kable(digits = 5)| iteration | lambda | log(likelihood) | score | hessian | diff(loglik) |
|---|---|---|---|---|---|
| 0 | 3.00000 | -40.0610 | 11.00000 | -8.00000 | NA |
| 1 | 4.37500 | -30.7708 | 3.45714 | -3.76163 | 9.29018 |
| 2 | 5.29405 | -28.9897 | 0.60016 | -2.56895 | 1.78110 |
| 3 | 5.52768 | -28.9176 | 0.02537 | -2.35639 | 0.07210 |
| 4 | 5.53844 | -28.9175 | 0.00005 | -2.34724 | 0.00014 |
| 5 | 5.53846 | -28.9175 | 0.00000 | -2.34722 | 0.00000 |
The final estimate matches the closed-form MLE, \(\bar x = 5.53846\) (Exercise 23).
Show R code
3 Maximum likelihood for univariate Gaussian models
Suppose \(X_1, \ldots, X_n \ \sim_{\operatorname{iid}}\ \operatorname{N}\mathopen{}\left(\mu, \sigma^2\right)\mathclose{}\), and let \(x_1, \ldots, x_n\) be the observed values. The parameter vector is \(\tilde{\theta}= (\mu, \sigma^2)\). We treat \(\sigma^2\), rather than \(\sigma\), as the second parameter, and differentiate with respect to \(\sigma^2\) directly.
By Theorem 1, the likelihood is:
\[ \mathcal{L}(\mu, \sigma^2) = \prod_{i=1}^n (2\pi\sigma^2)^{-1/2} \operatorname{exp}\mathopen{}\left\{-\frac{(x_i - \mu)^2}{2\sigma^2}\right\}\mathclose{} \]
and by Theorem 4, the log-likelihood is:
\[ \begin{aligned} \ell(\mu, \sigma^2) &= \sum_{i=1}^n \mathopen{}\left(-\frac{1}{2}\operatorname{log}\mathopen{}\left\{2\pi\sigma^2\right\}\mathclose{} - \frac{(x_i - \mu)^2}{2\sigma^2}\right)\mathclose{} && \text{(log of each Gaussian density)}\\ &= -\frac{n}{2}\operatorname{log}\mathopen{}\left\{2\pi\right\}\mathclose{} - \frac{n}{2}\operatorname{log}\mathopen{}\left\{\sigma^2\right\}\mathclose{} - \frac{1}{2\sigma^2}\sum_{i=1}^n (x_i - \mu)^2 && \text{(split the sum; log of a product)} \end{aligned} \]
3.1 The score function
The score function is the vector of the two partial derivatives:
\[ \ell'(\mu, \sigma^2) = \begin{pmatrix} \frac{\partial}{\partial \mu}\ell(\mu, \sigma^2) \\ \frac{\partial}{\partial \sigma^2}\ell(\mu, \sigma^2) \end{pmatrix} \]
For the first entry, only the last term of \(\ell\) depends on \(\mu\):
\[ \begin{aligned} \frac{\partial}{\partial \mu}\ell &= -\frac{1}{2\sigma^2}\sum_{i=1}^n \frac{\partial}{\partial \mu}(x_i - \mu)^2 && \text{(linearity of differentiation)}\\ &= -\frac{1}{2\sigma^2}\sum_{i=1}^n -2(x_i - \mu) && \text{(chain rule)}\\ &= \frac{1}{\sigma^2}\mathopen{}\left(\sum_{i=1}^n x_i - n\mu\right)\mathclose{} && \text{(simplify; split the sum)} \end{aligned} \]
For the second entry, write \(\sigma^2\) as a single variable \(v\), so that \(\ell\) contains \(-\frac{n}{2}\log v\) and \(-\frac{1}{2}v^{-1}\sum_i (x_i - \mu)^2\):
\[ \begin{aligned} \frac{\partial}{\partial \sigma^2}\ell &= -\frac{n}{2}\mathopen{}\left(\sigma^2\right)\mathclose{}^{-1} + \frac{1}{2}\mathopen{}\left(\sigma^2\right)\mathclose{}^{-2}\sum_{i=1}^n (x_i - \mu)^2 && \text{(derivatives of $\log v$ and $v^{-1}$)} \end{aligned} \]
3.2 MLE of \(\mu\)
Setting \(\frac{\partial}{\partial \mu}\ell = 0\):
\[ \begin{aligned} 0 &= \frac{1}{\sigma^2}\mathopen{}\left(\sum_{i=1}^n x_i - n\mu\right)\mathclose{} && \text{(score for $\mu$)}\\ n\mu &= \sum_{i=1}^n x_i && \text{(multiply by $\sigma^2$; add $n\mu$)}\\ \mu &= \bar x && \text{(divide by $n$)} \end{aligned} \]
This solution does not depend on \(\sigma^2\). The second derivative is
\[\frac{\partial^2 \ell}{\partial \mu^2} = \frac{\partial}{\partial \mu}\frac{1}{\sigma^2}\mathopen{}\left(\sum_{i=1}^n x_i - n\mu\right)\mathclose{} = -\frac{n}{\sigma^2} < 0,\]
so for every fixed \(\sigma^2\), \(\ell\) is maximized over \(\mu\) at \(\bar x\), and \(\hat\mu_{\text{ML}} = \bar x\).
3.3 MLE of \(\sigma^2\)
Setting \(\frac{\partial}{\partial \sigma^2}\ell = 0\):
\[ \begin{aligned} 0 &= -\frac{n}{2}\mathopen{}\left(\sigma^2\right)\mathclose{}^{-1} + \frac{1}{2}\mathopen{}\left(\sigma^2\right)\mathclose{}^{-2}\sum_{i=1}^n (x_i - \mu)^2 && \text{(score for $\sigma^2$)}\\ \frac{n}{2}\mathopen{}\left(\sigma^2\right)\mathclose{}^{-1} &= \frac{1}{2}\mathopen{}\left(\sigma^2\right)\mathclose{}^{-2}\sum_{i=1}^n (x_i - \mu)^2 && \text{(add $\tfrac{n}{2}(\sigma^2)^{-1}$)}\\ \sigma^2 &= \frac{1}{n}\sum_{i=1}^n (x_i - \mu)^2 && \text{(multiply by $2(\sigma^2)^2/n$)} \end{aligned} \]
Substituting the maximizer \(\mu = \bar x\), which does not depend on \(\sigma^2\), maximizes the profile log-likelihood of Example 13, and gives:
\[\hat{\sigma}^2_{\text{ML}} = \frac{1}{n}\sum_{i=1}^n (x_i - \bar x)^2\]
The profile log-likelihood, \(\ell_p(\sigma^2) = -\frac{n}{2}\operatorname{log}\mathopen{}\left\{2\pi\right\}\mathclose{} - \frac{n}{2}\operatorname{log}\mathopen{}\left\{\sigma^2\right\}\mathclose{} - \frac{n\hat\sigma^2_{\text{ML}}}{2\sigma^2}\), increases for \(\sigma^2 < \hat\sigma^2_{\text{ML}}\) and decreases for \(\sigma^2 > \hat\sigma^2_{\text{ML}}\), because its derivative, \(\frac{n}{2}\mathopen{}\left(\sigma^2\right)\mathclose{}^{-2}\mathopen{}\left(\hat\sigma^2_{\text{ML}} - \sigma^2\right)\mathclose{}\), has the sign of \(\hat\sigma^2_{\text{ML}} - \sigma^2\). So \((\bar x, \hat\sigma^2_{\text{ML}})\) is the global maximizer, provided the \(x_i\) are not all equal.
Differentiating with respect to \(\sigma^2\) as a single variable, rather than with respect to \(\sigma\), keeps the algebra short. Replacing \(\sigma^2\) with the precision \(\tau \stackrel{\text{def}}{=}1/\sigma^2\) and differentiating with respect to \(\tau\) can be shorter still, because \(\tau\) enters the log-likelihood as \(\frac{n}{2}\log\tau - \frac{\tau}{2}\sum_{i=1}^n (x_i - \mu)^2\). By the invariance of maximum likelihood estimates, \(\hat\tau_{\text{ML}} = 1/\hat\sigma^2_{\text{ML}}\).
This MLE divides by \(n\), so it is a biased estimator of \(\sigma^2\) (bias of the divide-by-\(n\) estimator).
3.4 Second derivatives
The remaining second derivatives are:
\[ \begin{aligned} \frac{\partial^2 \ell}{\partial (\sigma^2)^2} &= \frac{\partial}{\partial \sigma^2}\mathopen{}\left(-\frac{n}{2}\mathopen{}\left(\sigma^2\right)\mathclose{}^{-1} + \frac{1}{2}\mathopen{}\left(\sigma^2\right)\mathclose{}^{-2}\sum_{i=1}^n (x_i - \mu)^2\right)\mathclose{} && \text{(differentiate the score for $\sigma^2$)}\\ &= \frac{n}{2}\mathopen{}\left(\sigma^2\right)\mathclose{}^{-2} - \mathopen{}\left(\sigma^2\right)\mathclose{}^{-3}\sum_{i=1}^n (x_i - \mu)^2 && \text{(power rule)}\\ \frac{\partial^2 \ell}{\partial \mu \, \partial \sigma^2} &= \frac{\partial}{\partial \mu}\mathopen{}\left(\frac{1}{2}\mathopen{}\left(\sigma^2\right)\mathclose{}^{-2}\sum_{i=1}^n (x_i - \mu)^2\right)\mathclose{} && \text{(only this term depends on $\mu$)}\\ &= -\mathopen{}\left(\sigma^2\right)\mathclose{}^{-2}\sum_{i=1}^n (x_i - \mu) && \text{(chain rule)} \end{aligned} \]
At the MLE, \(\sum_{i=1}^n (x_i - \bar x) = 0\) and \(\sum_{i=1}^n (x_i - \bar x)^2 = n\hat\sigma^2\) (writing \(\hat\sigma^2\) for \(\hat\sigma^2_{\text{ML}}\)), so:
\[ \begin{aligned} \frac{\partial^2 \ell}{\partial (\sigma^2)^2}\bigg|_{\text{MLE}} &= \frac{n}{2}\mathopen{}\left(\hat\sigma^2\right)\mathclose{}^{-2} - \mathopen{}\left(\hat\sigma^2\right)\mathclose{}^{-3} n\hat\sigma^2 && \text{(substitute)}\\ &= -\frac{n}{2}\mathopen{}\left(\hat\sigma^2\right)\mathclose{}^{-2} && \text{(simplify)}\\ \frac{\partial^2 \ell}{\partial \mu \, \partial \sigma^2}\bigg|_{\text{MLE}} &= 0 && \text{(substitute)} \end{aligned} \]
3.5 Information matrix and standard errors
Collecting the second derivatives at the MLE, the observed information is
\[ I(\hat\mu, \hat\sigma^2) = \begin{bmatrix} \frac{n}{\hat\sigma^2} & 0 \\ 0 & \frac{n}{2\mathopen{}\left(\hat\sigma^2\right)\mathclose{}^2} \end{bmatrix} \]
Its diagonal entries are positive and its off-diagonal entries are zero, so it is positive definite, consistent with the MLE being a maximum. The inverse of a diagonal matrix inverts each diagonal entry, so:
\[ \mathopen{}\left(I(\hat\mu, \hat\sigma^2)\right)^{-1}\mathclose{} = \begin{bmatrix} \frac{\hat\sigma^2}{n} & 0 \\ 0 & \frac{2\mathopen{}\left(\hat\sigma^2\right)\mathclose{}^2}{n} \end{bmatrix} \]
By Theorem 9, the estimated standard errors are \(\mathop{\widehat{\operatorname{SE}}}\nolimits\mathopen{}\left(\hat\mu\right)\mathclose{} = \hat\sigma/\sqrt{n}\) and \(\mathop{\widehat{\operatorname{SE}}}\nolimits\mathopen{}\left(\hat\sigma^2\right)\mathclose{} = \hat\sigma^2 \sqrt{2/n}\), and the two estimates are approximately uncorrelated in large samples.
See also (Casella and Berger 2002, Example 7.2.12).
4 Example: hormone therapy study
This example fits a Gaussian model to real data by maximum likelihood, and then uses simulation to examine the properties of maximum likelihood estimation for that model.
4.1 Data
The “heart and estrogen/progestin study” (HERS) was a clinical trial of hormone therapy for prevention of recurrent heart attacks and death among 2,763 post-menopausal women with existing coronary heart disease (CHD) (Hulley et al. 1998).
The trial was conducted at 20 US clinical centers. Participants were randomized to receive either conjugated equine estrogens (0.625 mg/day) plus medroxyprogesterone acetate (2.5 mg/day) or a matching placebo (Hulley et al. 1998). Women were followed for an average of 4.1 years (Hulley et al. 1998).
The primary outcome was nonfatal myocardial infarction or CHD death (Hulley et al. 1998).
We model the distribution of fasting glucose among HERS participants who do not have diabetes and do not exercise.
The HERS data are distributed with Vittinghoff et al. (2012) on the book’s companion website:
# one unbroken string, so that link checkers test the whole URL:
url <- "https://regression.ucsf.edu/sites/g/files/tkssra16191/files/wysiwyg/home/data/hersdata.dta" # nolint: line_length_linter.
hers <- haven::read_dta(url)The rmb R package includes the same file, which these notes use so that rendering does not depend on the website (Table 6):
hers <- rmb::hers |> haven::zap_labels()Show R code
hers |> head()To keep the likelihood graphs readable, we use only the first 100 eligible participants (Figure 14); with the whole subset, the likelihood would be too concentrated to graph clearly.
Show R code
plot1 <-
data1 |>
ggplot2::ggplot() +
ggplot2::aes(x = glucose) +
ggplot2::geom_histogram(
ggplot2::aes(y = ggplot2::after_stat(density)),
bins = 20
) +
ggplot2::xlab("Fasting glucose (mg/dL)")
print(plot1)
The histogram is irregular, as histograms of 100 observations often are, with a somewhat longer right tail than left tail. A Gaussian model is a rough but usable starting point.
4.2 Maximum likelihood estimates
By the Gaussian MLEs, \(\hat\mu_{\text{ML}} = \bar x\) and \(\hat\sigma^2_{\text{ML}} = \frac{1}{n}\sum_i (x_i - \bar x)^2\):
Figure 15 superimposes the fitted Gaussian density on the histogram.
Show R code
plot1 +
ggplot2::geom_function(
fun = function(x) dnorm(x, mean = mu_hat, sd = sigma_hat),
col = "red"
)
The fitted curve follows the overall shape of the histogram, but it underestimates the frequency of values between 110 and 122 mg/dL, consistent with the histogram’s longer right tail.
4.3 Likelihood and log-likelihood functions
It is numerically better to compute the log-likelihood first and exponentiate it to get the likelihood, because a product of 100 densities can underflow to zero:
loglik <- function(mu, sigma, x) {
n <- length(x)
normalizing_constant <- -n / 2 * log(2 * pi * sigma^2)
# written with sum(x), sum(x^2) so that `mu` can be a vector of values:
kernel <- -1 / (2 * sigma^2) * (sum(x^2) - 2 * sum(x) * mu + n * mu^2)
normalizing_constant + kernel
}
lik <- function(...) exp(loglik(...))Figure 16 graphs the likelihood and log-likelihood as functions of \(\mu\), with \(\sigma\) fixed at \(\hat\sigma_{\text{ML}}\).
Show R code
ggplot2::ggplot() +
ggplot2::geom_function(
fun = lik,
args = list(sigma = sigma_hat, x = glucose_data)
) +
ggplot2::xlim(mu_hat + c(-1, 1) * sigma_hat) +
ggplot2::xlab("mu") +
ggplot2::ylab("likelihood") +
ggplot2::geom_vline(xintercept = mu_hat, col = "red")
Show R code
ggplot2::ggplot() +
ggplot2::geom_function(
fun = loglik,
args = list(sigma = sigma_hat, x = glucose_data)
) +
ggplot2::xlim(mu_hat + c(-1, 1) * sigma_hat) +
ggplot2::xlab("mu") +
ggplot2::ylab("log-likelihood") +
ggplot2::geom_vline(xintercept = mu_hat, col = "red")
Figure 17 graphs them as functions of \(\sigma\), with \(\mu\) fixed at \(\hat\mu_{\text{ML}}\).
Show R code
ggplot2::ggplot() +
ggplot2::geom_function(
fun = lik,
args = list(mu = mu_hat, x = glucose_data)
) +
ggplot2::xlim(sigma_hat * c(0.9, 1.1)) +
ggplot2::geom_vline(xintercept = sigma_hat, col = "red") +
ggplot2::xlab("sigma") +
ggplot2::ylab("likelihood")
Show R code
ggplot2::ggplot() +
ggplot2::geom_function(
fun = loglik,
args = list(mu = mu_hat, x = glucose_data)
) +
ggplot2::xlim(sigma_hat * c(0.9, 1.1)) +
ggplot2::geom_vline(xintercept = sigma_hat, col = "red") +
ggplot2::xlab("sigma") +
ggplot2::ylab("log-likelihood")
4.4 Log-likelihood surface
Figure 18 graphs the log-likelihood over both parameters at once.
Show R code
n_points <- 25
mu_grid <- seq(94, 104, length.out = n_points)
sigma_grid <- seq(7, 15, length.out = n_points)
lliks <- outer(mu_grid, sigma_grid, loglik, x = glucose_data)
plotly::plot_ly(
type = "surface",
x = ~mu_grid,
y = ~sigma_grid,
# plot_ly() expects rows of z to correspond to y values,
# so the matrix is transposed;
# see https://stackoverflow.com/questions/69472185
z = ~ t(lliks)
) |>
plotly::layout(
scene = list(
xaxis = list(title = "mu"),
yaxis = list(title = "sigma"),
zaxis = list(title = "log-likelihood")
)
)4.5 Standard errors by sample size
By Section 3.5, the estimated standard error of \(\hat\mu_{\text{ML}}\) is
\[ \mathop{\widehat{\operatorname{SE}}}\nolimits\mathopen{}\left(\hat\mu\right)\mathclose{} = \sqrt{\mathopen{}\left[\mathopen{}\left(I(\hat\mu, \hat\sigma^2)\right)^{-1}\mathclose{}\right]\mathclose{}_{11}} = \frac{\hat\sigma}{\sqrt{n}} \]
which shrinks in proportion to \(1/\sqrt{n}\) (Figure 19).
Show R code
se_mu_hat <- function(n, sigma) sigma / sqrt(n)
ggplot2::ggplot() +
ggplot2::geom_function(fun = se_mu_hat, args = list(sigma = sigma_hat)) +
ggplot2::scale_x_log10(
limits = c(10, 10^5), name = "Sample size",
labels = scales::label_comma()
) +
ggplot2::ylab("Standard error of mu-hat (mg/dL)")
4.6 Power
Suppose we test the null hypothesis \(H_0: \mu = \mu_0\), with \(\mu_0 = 95\) mg/dL, at significance level \(\alpha = 0.05\), and suppose for simplicity that \(\sigma\) is known, equal to \(\hat\sigma_{\text{ML}}\). Then under \(H_0\), \(\bar X \sim \operatorname{N}\mathopen{}\left(\mu_0, \sigma^2/n\right)\mathclose{}\), and the test rejects \(H_0\) when \(\bar x\) falls outside the non-rejection interval
\[\mu_0 \pm z_{1 - \alpha/2} \frac{\sigma}{\sqrt{n}}\]
For this test, under \(\mu = \mu_1\), \(\bar X \sim \operatorname{N}\mathopen{}\left(\mu_1, \sigma^2/n\right)\mathclose{}\), so:
\[ \text{power}(\mu_1) = \Phi\mathopen{}\left(\frac{\mu_0 - z_{1-\alpha/2}\,\sigma/\sqrt{n} - \mu_1}{\sigma/\sqrt{n}}\right)\mathclose{} + 1 - \Phi\mathopen{}\left(\frac{\mu_0 + z_{1-\alpha/2}\,\sigma/\sqrt{n} - \mu_1}{\sigma/\sqrt{n}}\right)\mathclose{} \]
where \(\Phi\) is the standard Gaussian CDF. For example, the power against \(\mu_1 = 100\) mg/dL is:
Figure 20 graphs the power as a function of the sample size.
Show R code
The alternative \(\mu_1\) should be chosen before seeing the data, as a difference worth detecting. Power computed at \(\mu_1 = \hat\mu\) (“observed power”) is a function of the p-value, so it adds no information about the data already analyzed (Hoenig and Heisey 2001).
4.7 Simulation
To check how maximum likelihood estimation behaves for this model, we simulate many datasets from a Gaussian distribution whose parameters equal the HERS estimates, analyze each one, and summarize the results.
do_one_sim() simulates and analyzes one dataset: it computes \(\hat\mu\), its estimated standard error, a 95% \(t\)-based confidence interval for \(\mu\), and the \(t\)-test of \(H_0: \mu = \mu_0\).
do_one_sim <- function(n, mu, mu0, sigma2, return_data = FALSE) {
# generate data
x <- rnorm(n = n, mean = mu, sd = sqrt(sigma2))
# analyze data
est <- mean(x)
se_est <- sd(x) / sqrt(n)
confint <- est + c(-1, 1) * se_est * qt(0.975, df = n - 1)
tstat <- abs(est - mu0) / se_est
pval <- 2 * pt(q = tstat, df = n - 1, lower.tail = FALSE)
results <- tibble::tibble(
mu_hat = est,
sigma_hat = sd(x),
se_hat = se_est,
confint_left = confint[1],
confint_right = confint[2],
tstat = tstat,
pval = pval,
confint_covers = dplyr::between(mu, confint[1], confint[2]),
test_rejects = pval < 0.05
)
if (return_data) {
list(data = x, results = results)
} else {
results
}
}To check do_one_sim(), we compare its output with stats::t.test() on the same simulated data:
set.seed(1)
sim_output <- do_one_sim(
n = 100, mu = mu_hat, mu0 = 80, sigma2 = sigma_sq_hat,
return_data = TRUE
)
t_test <- t.test(sim_output$data, mu = 80)
dplyr::bind_rows(
`do_one_sim()` = sim_output$results |>
dplyr::select(mu_hat, se_hat, confint_left, confint_right, pval),
`t.test()` = tibble::tibble(
mu_hat = unname(t_test$estimate),
se_hat = t_test$stderr,
confint_left = t_test$conf.int[1],
confint_right = t_test$conf.int[2],
pval = t_test$p.value
),
.id = "source"
)The two rows agree.
do_n_sims() repeats the simulation n_sims times, with a different random seed for each dataset:
do_n_sims <- function(n_sims = 1000, ...) {
lapply(seq_len(n_sims), function(i) {
set.seed(i)
do_one_sim(...) |>
dplyr::mutate(sim_number = i, .before = dplyr::everything())
}) |>
dplyr::bind_rows()
}
sim_results <- do_n_sims(
n_sims = 1000,
n = 100, mu = mu_hat, mu0 = 0.9 * mu_hat, sigma2 = sigma_sq_hat
)
sim_resultssummarize_sim() compares the simulation results with the true data-generating parameters:
summarize_sim <- function(sim_results, mu, sigma2, n) {
true_se <- sqrt(sigma2 / n)
tibble::tibble(
bias_mu_hat = mean(sim_results$mu_hat) - mu,
sd_mu_hat = sd(sim_results$mu_hat),
true_se = true_se,
bias_se_hat = mean(sim_results$se_hat) - true_se,
coverage = mean(sim_results$confint_covers),
power = mean(sim_results$test_rejects)
)
}
sim_summary <- summarize_sim(
sim_results,
mu = mu_hat, sigma2 = sigma_sq_hat, n = 100
)
sim_summaryAcross 1000 simulated datasets:
- the average error of \(\hat\mu\) is -0.005 mg/dL, small relative to its standard error of 1.02 mg/dL, consistent with \(\hat\mu\) being unbiased;
- the standard deviation of the \(\hat\mu\) values, 0.996, is close to the true standard error;
- the 95% confidence intervals covered the true \(\mu\) in 95.9% of datasets, close to their nominal 95%;
- the test of \(H_0: \mu = 0.9\,\hat\mu\) rejected in 100% of datasets: with \(n = 100\), a 10% difference in the mean is easy to detect.
Changing the sample size, the true \(\mu\), or \(\sigma^2\) in do_n_sims() shows how these properties depend on them.
5 Practice exercises
References
Casella, George, and Roger Berger. 2002. Statistical Inference. 2nd ed. Cengage Learning. https://www.cengage.com/c/statistical-inference-2e-casella-berger/9780534243128/.
Dobson, Annette J, and Adrian G Barnett. 2018. An Introduction to Generalized Linear Models. 4th ed. CRC press. https://doi.org/10.1201/9781315182780.
Dunn, Peter K, and Gordon K Smyth. 2018. Generalized Linear Models with Examples in R. Vol. 53. Springer. https://doi.org/10.1007/978-1-4419-0118-7.
Efron, Bradley, and David V Hinkley. 1978. “Assessing the Accuracy of the Maximum Likelihood Estimator: Observed Versus Expected Fisher Information.” Biometrika 65 (3): 457–83. https://doi.org/10.1093/biomet/65.3.457.
Hoenig, John M., and Dennis M. Heisey. 2001. “The Abuse of Power: The Pervasive Fallacy of Power Calculations for Data Analysis.” The American Statistician 55 (1): 19–24. https://doi.org/10.1198/000313001300339897.
Hogg, Robert V., Elliot A. Tanis, and Dale L. Zimmerman. 2019. Probability and Statistical Inference. Tenth edition. Pearson.
Hulley, Stephen, Deborah Grady, Trudy Bush, et al. 1998. “Randomized Trial of Estrogen Plus Progestin for Secondary Prevention of Coronary Heart Disease in Postmenopausal Women.” JAMA : The Journal of the American Medical Association (Chicago, IL) 280 (7): 605–13. https://doi.org/10.1001/jama.280.7.605.
Kleinbaum, David G, Lawrence L Kupper, Azhar Nizam, K Muller, and ES Rosenberg. 2014. Applied Regression Analysis and Other Multivariable Methods. 5th ed. Cengage Learning. https://www.cengage.com/c/applied-regression-analysis-and-other-multivariable-methods-5e-kleinbaum/9781285051086/.
Lehmann, E. L. 1999. Elements of Large-Sample Theory. Springer Texts in Statistics. Springer. https://doi.org/10.1007/b98855.
McLachlan, Geoffrey J, and Thriyambakam Krishnan. 2007. The EM Algorithm and Extensions. 2nd ed. John Wiley & Sons. https://doi.org/10.1002/9780470191613.
Newey, Whitney K, and Daniel McFadden. 1994. “Large Sample Estimation and Hypothesis Testing.” In Handbook of Econometrics, edited by Robert Engle and Dan McFadden, vol. 4. Elsevier. https://doi.org/10.1016/S1573-4412(05)80005-4.
Vittinghoff, Eric, David V Glidden, Stephen C Shiboski, and Charles E McCulloch. 2012. Regression Methods in Biostatistics: Linear, Logistic, Survival, and Repeated Measures Models. 2nd ed. Springer. https://doi.org/10.1007/978-1-4614-1353-0.
Wilks, Samuel S. 1938. “The Large-Sample Distribution of the Likelihood Ratio for Testing Composite Hypotheses.” The Annals of Mathematical Statistics 9 (1): 60–62. https://doi.org/10.1214/aoms/1177732360.
Wood, Simon N. 2017. Generalized Additive Models: An Introduction with r. 2nd ed. Chapman; Hall/CRC. https://doi.org/10.1201/9781315370279.