Maximum Likelihood Estimation

Maximum Likelihood Estimation in Health Economic Models

Maximum likelihood estimation picks, from all the parameter values an assumed probability model allows, the ones under which the data actually collected would have been most likely to arise. In health economics it is the routine way of fitting the parametric survival curves that are extrapolated beyond trial follow-up and the regression equations that turn patient characteristics into event risks. The same calculation returns standard errors and a covariance matrix, so a likelihood fit can feed probabilistic sensitivity analysis directly. R. A. Fisher gave the method its name in a paper published in April 1922, a decade after a first version in his 1912 paper on fitting frequency curves.

Likelihood as a function of the parameters

Fisher described the method as choosing the set of parameter values that makes the probability of the observed sample a maximum. Read as a function of the unknown parameters with the data held fixed, that probability is the likelihood. For independent observations it is a product, and in practice the logarithm is maximised because a sum is easier to handle and has the same maximum.

$$L(\theta; y) = \prod_{i=1}^{n} f_i(y_i; \theta), \qquad \log L(\theta; y) = \sum_{i=1}^{n} \log f_i(y_i; \theta)$$

where $\theta$ is the vector of unknown parameters, $y = (y_1, \dots, y_n)$ are the $n$ observed values, $f_i(y_i; \theta)$ is the probability or density of observation $i$ under the model, $L$ is the likelihood and $\log L$ is the log-likelihood. The maximum likelihood estimate $\hat\theta$ is the value with $\log L(\hat\theta; y) \ge \log L(\theta; y)$ for every allowed $\theta$.

Fisher was explicit that likelihood is not used as a synonym for probability. It expresses the relative frequencies with which each candidate value would yield the observed sample, and the likelihoods of all candidate values do not sum to one. He also stressed that the position of the maximum is entirely unchanged when the parameter is put on a different scale. This invariance lets an analyst estimate a hazard rate and convert it to a cycle probability without refitting.

Finding the maximum with the score and the information

The first derivative of the log-likelihood is the score. When the log-likelihood is concave, setting the score to zero and solving gives the estimate, and the curvature at the peak measures how precisely the data pin it down.

$$u(\theta) = \frac{\partial \log L(\theta; y)}{\partial \theta}, \qquad u(\hat\theta) = 0, \qquad \widehat{\text{var}}(\hat\theta) = I(\hat\theta)^{-1}$$

where $u(\theta)$ is the score vector, $I(\theta)$ is the information matrix, either minus the expected second derivatives of the log-likelihood or, as the observed information, minus the second derivatives evaluated with the data in hand, and $\widehat{\text{var}}(\hat\theta)$ is the estimated variance-covariance matrix of the estimates. A sharply peaked log-likelihood gives small standard errors and a flat one gives large ones.

Under regularity conditions the estimator is approximately normal in large samples, with mean equal to the true value and variance given by the inverse of the information. That result is the basis for the usual standard error, the Wald test and the large-sample confidence interval, covered in more depth on the page on asymptotic normality. The conditions include a true value inside the parameter space, not on its edge. Where the score equations have no closed-form solution, as for a Weibull curve or a logistic regression, software iterates. The Newton-Raphson method updates a trial value using the score and the matrix of second derivatives, and Fisher scoring replaces that matrix with its expected value.

Censored survival data and the likelihood

Trial time-to-event data are a common setting for maximum likelihood in health economics, because many patients are still alive or event-free when follow-up ends. A patient whose event is observed at time $t_i$ contributes the density at $t_i$. A patient who is censored at $t_i$ contributes only the probability of surviving beyond $t_i$, which is valid if the censoring itself carries no information about that patient's prospects.

$$\log L = \sum_{i=1}^{n} \left[ d_i \log h(t_i) - H(t_i) \right]$$

where $t_i$ is the follow-up time of patient $i$, $d_i$ is 1 if the event was observed and 0 if the time is censored, $h(t)$ is the hazard function and $H(t)$ is the cumulative hazard, so that survival is $S(t) = e^{-H(t)}$.

NICE Decision Support Unit Technical Support Document 14 (TSD 14) says that exponential, Weibull, Gompertz, log-logistic, log normal and generalised gamma models should all be considered when extrapolating trial survival with patient-level data. Each candidate is fitted by maximising this likelihood, and the maximised values are then compared through the Akaike information criterion and the Bayesian information criterion, alongside visual fit and clinical plausibility. The page on the parametric survival model covers the curve shapes and model choice.

Worked example with a constant hazard from trial follow-up

All numbers in this example are illustrative. A trial arm of 200 patients records 50 deaths, and total follow-up across all 200, censored patients included, is 400 person-years. The analyst fits an exponential model, which assumes a constant hazard rate $\lambda$, so $h(t) = \lambda$ and $H(t) = \lambda t$.

1. Write the log-likelihood. With $D$ deaths and total time at risk $T$, the censored-data log-likelihood reduces to the following.

$$\log L(\lambda) = D \log \lambda - \lambda T = 50 \log \lambda - 400 \lambda$$

where $D = 50$ is the number of deaths, $T = 400$ is the total person-years and $\lambda$ is the hazard per person-year.

2. Solve the score equation. The score is $D/\lambda - T$, and setting it to zero gives $\hat\lambda = D/T = 50/400 = 0.125$ deaths per person-year, the familiar deaths divided by exposure time.

3. Check that it is the peak. At $\lambda = 0.10$ the log-likelihood is 50 × (minus 2.3026) minus 40 = minus 155.13, at 0.125 it is 50 × (minus 2.0794) minus 50 = minus 153.97, and at 0.15 it is 50 × (minus 1.8971) minus 60 = minus 154.86, so 0.125 is highest.

4. Standard error from the information. The observed information is $D/\lambda^2$, which at the estimate is 50 / 0.015625 = 3,200. The variance is 1 / 3,200 = 0.0003125, so the standard error is 0.0177, the same as $\sqrt{50}/400$.

5. Interval on the log scale. A hazard cannot be negative, so the interval is better built for $\log \lambda$. Writing $\eta = \log \lambda$, the log-likelihood is $D\eta - T e^{\eta}$ and its second derivative at the estimate is minus 400 × 0.125 = minus 50, so the standard error of $\hat\eta$ is one divided by $\sqrt{50}$, or 0.1414. The interval is $\log 0.125$ = minus 2.0794 plus or minus 1.96 × 0.1414 = 0.2772, which gives minus 2.3566 to minus 1.8023 and, after exponentiating, a hazard of 0.0947 to 0.1649 per person-year.

6. Convert for a yearly cycle. By invariance the estimate of the annual transition probability is one minus $e^{-0.125}$, or 0.1175, and converting the interval limits the same way gives 0.0904 to 0.1520. The crude 50/200 = 0.25 ignores how long each patient was followed.

7. Interpretation. A Markov cohort model with a one-year cycle would use 0.1175 as its base-case annual death probability, and for probabilistic analysis it could draw $\log \lambda$ from a normal distribution with mean minus 2.0794 and standard deviation 0.1414. The constant-hazard assumption implies a mean survival of 1 / 0.125 = 8 years, which depends wholly on extrapolating the exponential shape well beyond an average follow-up of 2 years per patient (400 person-years across 200 patients).

Fitted parameters as decision model inputs

The ISPOR-SMDM task force report on parameter estimation notes that estimation returns a point estimate, a measure of precision and, in a multivariable fit, the covariance between estimates. For regression-based parameters it states that the covariance matrix defines the uncertainty in the coefficients, that multivariate normality is appropriate for the linear predictor, and that this can set the range for deterministic sensitivity analysis or form the basis for probabilistic sensitivity analysis. TSD 14 says the same for survival curves fitted to patient-level data, using the variance-covariance matrix of each fitted model. Correlated draws can be generated with a Cholesky decomposition of that matrix, so that the shape and scale of a Weibull curve move together as the fit implies.

TSD 14 also treats the choice between fitted distributions as a form of structural sensitivity analysis, which sampling from one covariance matrix does not capture. Guyot and colleagues note that individual patient data allow a range of parametric models to be investigated, and for trials that publish only Kaplan-Meier curves they set out an algorithm that derives a close approximation to those data. They judged it reliable for cost-effectiveness analysis, although hazard ratios were reasonably accurate only when numbers at risk or total events were reported.

The same machinery fits the risk equations in many models. A logistic regression for an event probability, a generalised linear model for costs and a parametric survival model with a treatment covariate are all estimated by maximising a likelihood, and their reported standard errors come from its curvature.

Assumptions behind a likelihood estimate

The method inherits every assumption written into the probability model. In the exponential example the estimator $D/T$ is optimal in the maximum likelihood sense only if the risk really is constant over time; if the hazard rises with age or time since diagnosis, the estimate averages over a pattern the model cannot represent. Observations are treated as independent given the model, and censoring is assumed to tell nothing about a patient's later prospects.

The standard errors rest on large-sample approximations. With few events, or a parameter that the data barely identify, the log-likelihood can be lopsided or flat and a normal-based interval misleading; a transformed scale, as for $\log \lambda$ above, helps when the parameter has a natural bound. Iterative fitting can also stop at a local peak or fail to converge, so starting values and convergence checks belong in the report.

Maximum likelihood beside moments and Bayesian methods

Fisher's 1922 paper set maximum likelihood against the method of moments, one of the two estimation methods that Fisher had found wanting in 1912, and much of its second half examines how efficient moment estimators are for families of curves. David Hand's commentary notes that Fisher had earlier objected that the method of moments chooses which moments to use without justification.

A Bayesian analysis multiplies the likelihood by a prior distribution and reports a posterior distribution. Gelman and colleagues note that the maximum likelihood estimate is the posterior mode under a uniform prior density, and that in large samples the prior becomes nearly irrelevant, so the two often agree. They part company when data are sparse or an informative prior is used, and only the Bayesian route gives direct probability statements about the parameter.

Misreadings and limits

Most errors with maximum likelihood come from treating the fitted model as more certain than it is. The estimate answers a narrow question about one assumed model and one dataset.

1. Reading the likelihood as the probability of the parameter. It ranks candidate values by how well they explain the data; a probability statement about the parameter needs a prior.
2. Comparing raw log-likelihoods across different distributions. TSD 14 records an appraisal in which model choice relied on the log-likelihood of non-nested curves and calls that approach flawed. The likelihood ratio test applies to nested models; non-nested curves are compared with information criteria.
3. Equating best fit with best extrapolation. TSD 14 stresses that the plausibility of the extrapolated portion matters even more than fit to the observed data, because it drives mean survival.
4. Feeding the point estimate alone into the model. Ignoring the covariance matrix leaves parameter uncertainty out of the analysis, and the ISPOR-SMDM report notes the argument that nonlinear models need probabilistic analysis to give appropriate expected values.

Sources

• Briggs AH, Weinstein MC, Fenwick EAL, Karnon J, Sculpher MJ, Paltiel AD. Model parameter estimation and uncertainty analysis: a report of the ISPOR-SMDM Modeling Good Research Practices Task Force Working Group-6. Medical Decision Making. 2012;32(5):722-732. https://doi.org/10.1177/0272989X12458348
• Fisher RA. On the mathematical foundations of theoretical statistics. Philosophical Transactions of the Royal Society of London, Series A. 1922;222:309-368. https://doi.org/10.1098/rsta.1922.0009
• Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, Rubin DB. Bayesian Data Analysis. 3rd ed. CRC Press; 2013. Chapter 4. http://www.stat.columbia.edu/~gelman/book/BDA3.pdf
• Guyot P, Ades AE, Ouwens MJNM, Welton NJ. Enhanced secondary analysis of survival data: reconstructing the data from published Kaplan-Meier survival curves. BMC Medical Research Methodology. 2012;12:9. https://doi.org/10.1186/1471-2288-12-9
• Hand DJ. From evidence to understanding: a commentary on Fisher (1922) 'On the mathematical foundations of theoretical statistics'. Philosophical Transactions of the Royal Society A. 2015;373(2039):20140252. https://doi.org/10.1098/rsta.2014.0252
• Latimer N. NICE DSU Technical Support Document 14: Survival analysis for economic evaluations alongside clinical trials, extrapolation with patient-level data. Sheffield: Decision Support Unit, ScHARR; 2011, last updated March 2013. https://doi.org/10.15131/shef.data.33745843
• Rodríguez G. Generalized Linear Models, lecture notes. Appendix A, Review of likelihood theory (revised November 2001), and Chapter 7, Survival models (revised September 2010). https://grodri.github.io/glms/notes/a1.pdf and https://grodri.github.io/glms/notes/c7.pdf