Constructing a Profile Likelihood

Why Reoptimize at Each Fixed Value?

A fitted QSP model can reproduce observed data through compensation among parameters: changing one parameter may be offset by changing another. To examine the likelihood associated with a particular value of a parameter, hold that value fixed and let the other parameters adjust as well as they can. Repeating this operation over a grid constructs a profile likelihood. It is a computational procedure for comparing fitted likelihoods under constraints, not a new model or a substitute for specifying the observation model.

The calculation depends on the full likelihood and its parameter constraints. It also depends on actually finding the best fit of the remaining parameters at each grid point. A local optimization that stops at a poor solution gives a poor profile value, even if the optimization software reports convergence.

Definition and Calculation

Definition: Write the model parameter vector as $\theta=(\psi,\lambda)$, where $\psi$ is the scalar parameter of interest and $\lambda$ contains all other parameters. For likelihood $L(\psi,\lambda)$ and allowed nuisance-parameter set $\Lambda(\psi)$, the profile likelihood is $$ L_p(\psi)=\max_{\lambda\in\Lambda(\psi)} L(\psi,\lambda). $$ The profile log-likelihood is $\ell_p(\psi)=\log L_p(\psi)$. At each chosen value of $\psi$, the maximization is over every remaining free parameter, subject to the model's constraints.

The profile is defined relative to a specified likelihood: that likelihood encodes the model predictions, observation-error distribution, and any parameter-dependent terms. For example, if observations $y_i$ are conditionally independent and normally distributed with model means $f_i(\theta)$ and known standard deviations $\sigma_i$, then, up to a constant independent of $\theta$,

\ell(\theta)=-\frac{1}{2}\sum_{i=1}^{n} \left(\frac{y_i-f_i(\theta)}{\sigma_i}\right)^2.

This expression is appropriate only when that normal observation model and its independence assumptions are justified. If the error distribution or its parameters differ, use the corresponding likelihood instead. Maximizing the likelihood is equivalent to minimizing the weighted sum of squared residuals in this normal model.

A useful justification follows directly from constrained optimization. The full maximum likelihood estimate maximizes $L(\psi,\lambda)$ over both $\psi$ and $\lambda$. If the parameter space is the same and each constrained maximum is found, then maximizing the profile over $\psi$ recovers the full maximum:

\max_\psi L_p(\psi) =\max_\psi\max_{\lambda\in\Lambda(\psi)}L(\psi,\lambda) =\max_{(\psi,\lambda)}L(\psi,\lambda).

This identity is also a practical check: the largest profile value should agree with the full-fit likelihood, up to numerical tolerance. It can fail in a computed profile if a constrained optimization has not reached the relevant maximum, if constraints were implemented inconsistently, or if the grid missed the full-fit value.

Worked Example: Profiling the Half-Maximal Concentration in an Emax Model

Worked Example: Profiling K in a Gaussian Emax Model

Suppose a hypothetical experiment measures responses $y=(2,5)$ at concentrations $C=(1,4)$. Use the model $f(C)=E_{\max}C/(K+C)$, with independent normal errors of known standard deviation $\sigma=1$. Both parameters are positive. To construct profile points for $K$, fix a value of $K$ and reoptimize $E_{\max}$.

At fixed $K$, let $q_i=C_i/(K+C_i)$. The residual sum of squares is $\sum_i(y_i-E_{\max}q_i)^2$. Differentiating with respect to $E_{\max}$, setting the derivative to zero, and solving gives

\widehat E_{\max}(K)=\frac{\sum_i q_i y_i}{\sum_i q_i^2}.

For $K=0.5$, $q=(2/3,8/9)$. Therefore $$ \widehat E_{\max}(0.5)= \frac{(2/3)(2)+(8/9)(5)}{(2/3)^2+(8/9)^2} =\frac{52/9}{100/81}=4.68. $$ The fitted responses are $4.68(2/3)=3.12$ and $4.68(8/9)=4.16$. Thus $$ SSE_p(0.5)=(2-3.12)^2+(5-4.16)^2=1.2544+0.7056=1.96. $$

For comparison, at $K=1$, $q=(1/2,4/5)$, giving $$ \widehat E_{\max}(1)= \frac{(1/2)(2)+(4/5)(5)}{(1/2)^2+(4/5)^2} =\frac{5}{0.89}\approx5.618. $$ The fitted responses are approximately $2.809$ and $4.494$, so $SSE_p(1)\approx0.9101$. Under the stated normal model, the log-likelihood difference between these two profile points is $$ \ell_p(0.5)-\ell_p(1) =-\frac{1.96-0.9101}{2(1)^2} \approx-0.5249. $$ This calculation shows how to obtain and compare profile values; it does not by itself establish why the profile has this pattern or provide a confidence interval.

Worked Example: A Profile Point for a Nonlinear Decay Model

Worked Example: Fixing a Decay Rate and Reoptimizing Amplitude

Consider a hypothetical concentration-time model $f(t)=A\exp(-kt)$, with observations $y=(6.2,3.8,2.5)$ at times $t=(1,2,3)$ hours. Assume independent normal errors with known $\sigma=0.2$ concentration units and positive $A$ and $k$. Calculate the profile point at $k=0.5$ per hour.

With $q_i=\exp(-kt_i)$, the model is linear in the nuisance parameter $A$ at fixed $k$. The least-squares estimate is $\widehat A(k)=\sum_i q_i y_i/\sum_i q_i^2$. At $k=0.5$, $q\approx(0.60653,0.36788,0.22313)$, so $$ \widehat A(0.5)= \frac{(0.60653)(6.2)+(0.36788)(3.8)+(0.22313)(2.5)} {0.60653^2+0.36788^2+0.22313^2} \approx\frac{5.7163}{0.5530}=10.34. $$

The fitted concentrations are approximately $(6.270,3.803,2.306)$. The residual sum of squares is $$ SSE_p(0.5)\approx(6.2-6.270)^2+(3.8-3.803)^2+(2.5-2.306)^2 \approx0.0425. $$ Thus, up to a parameter-independent constant, the profile log-likelihood at this rate is $-0.0425/[2(0.2)^2]\approx-0.531$. Repeating the calculation at other values of $k$ supplies the grid of profile values. This example has an analytic nuisance-parameter solution; for a general QSP model, the remaining parameters may require numerical optimization.

Worked Example: Reoptimizing Multiple Nuisance Parameters

Worked Example: Baseline and Amplitude at a Fixed Rate

Suppose a hypothetical response follows $f(t)=b+A[1-\exp(-kt)]$, with observations $y=(1,3,4)$ at $t=(0,1,2)$. At a fixed $k$, both baseline $b$ and amplitude $A$ must be reoptimized. Assume independent normal errors with known standard deviation $\sigma=0.2$, and allow $b$ and $A$ to be real in this illustrative calculation.

At $k=\log(2)$ per time unit, the predictor values $q_i=1-\exp(-kt_i)$ are $(0,0.5,0.75)$. The data satisfy $y_i=1+4q_i$ exactly, so least squares gives $\widehat b=1$, $\widehat A=4$, fitted responses $(1,3,4)$, and $SSE_p(\log 2)=0$.

At $k=0.4$ per time unit, $q\approx(0,0.32968,0.55067)$. Regressing $y_i$ on an intercept and $q_i$ gives $$ \widehat A= \frac{\sum_i(q_i-\bar q)(y_i-\bar y)} {\sum_i(q_i-\bar q)^2} \approx\frac{0.8441}{0.1536}=5.496, \qquad \widehat b=\bar y-\widehat A\bar q\approx1.054. $$ The fitted responses are approximately $(1.054,2.866,4.080)$, giving $SSE_p(0.4)\approx0.0273$. This is a second profile point obtained by reoptimizing both nuisance parameters, rather than holding their estimates at the $k=\log(2)$ fit.

A Stable Computational Workflow

In a nonlinear dynamical model, each fixed grid value may require solving the ODE repeatedly during optimization. Numerical integration error, parameter scaling, constraints, and local optima can all affect the result. A disciplined workflow makes the calculation more reproducible and helps distinguish genuine model behavior from failed optimization.

Specify the likelihood and parameter domain.
State the observation model, fixed data and design, free parameters, and all bounds or other constraints. Keep the same likelihood and conventions at every grid point.

Obtain and check the full fit.
Fit all parameters jointly. Record the estimates, likelihood, optimizer status, and relevant numerical settings. The full-fit result provides a reference for the profile computation.

Choose a grid for the fixed parameter.
Cover the scientifically relevant allowed range. Use a log-spaced grid for a positive parameter spanning several orders of magnitude, or a linear grid when equal absolute increments are appropriate. Include the full-fit estimate when possible.

Optimize all remaining parameters at each grid point.
Fix the grid value exactly, then maximize the likelihood over every other free parameter. If using a minimizer, minimize the negative log-likelihood. Enforce the original bounds and evaluate the same observation model.

Check and refine the computations.
Save the objective, nuisance estimates, convergence information, and residual or prediction diagnostics at each point. Retry suspicious points with multiple starts or tighter numerical tolerances, and refine the grid where neighboring points or optimizer results require closer examination.

Parameter scaling matters because an optimizer's step sizes and stopping rules operate numerically, not in scientific units. If one rate is near $10^{-4}$ per hour and another parameter is near $10^2$, optimizing on their raw scales can make movement in one direction dominate. Rescaling parameters to comparable numerical ranges, or using log parameters for strictly positive quantities, can improve conditioning. When optimizing a log parameter, transform bounds and starting values consistently; the likelihood must still be evaluated using the corresponding original-scale parameter.

A warm start uses the nuisance estimates from a neighboring grid point as an initial guess. It often reduces computation when the optimum changes smoothly, but it can follow a local branch and miss a better solution. At important or suspicious points, supplement warm starts with several dispersed starting values, including starts near the full-fit estimates. Retain the highest likelihood among valid converged fits, not simply the first returned solution.

Diagnostics, Alternatives, and Limitations

A profile is a set of constrained optimization results, so inspect the optimization as well as the plotted values. Useful warning signs include abrupt changes in nuisance estimates accompanied by failed convergence, a point with a lower log-likelihood than a nearby point for no apparent numerical reason, large gradients at the reported solution, or parameters repeatedly landing on unintended bounds. Such signs do not prove an error, but they warrant checking the fit, tolerances, ODE solver, bounds, and initial values. Recompute questionable points using independent starts and, where appropriate, a finer grid.

Compare the best profile value with the unrestricted maximum likelihood. If the grid includes the full-fit estimate, its constrained fit should reproduce the unrestricted fit to numerical tolerance. If not, investigate whether the nuisance optimization failed or the constrained and full fits used different parameter domains, likelihood calculations, or numerical settings. A profile evaluated only on a coarse grid can also miss important locations between grid points, so grid refinement is part of computation rather than an optional cosmetic step.

Approach	What is optimized?	How it differs from profiling
Profile likelihood	At each fixed target value, all remaining free parameters	Accounts for nuisance-parameter adjustment under that constraint
Plug-in likelihood slice	No nuisance parameters; they remain fixed at chosen estimates	Can understate compensation because the nuisance parameters cannot adjust
Joint parameter grid	Evaluates combinations on a multidimensional grid	Can be costly as dimension grows; profiling uses optimization over nuisance dimensions
Wald approximation	Uses local curvature near the full fit	Does not perform a constrained refit at each target value

Profiling is more computationally demanding than evaluating a plug-in slice because it requires repeated optimization. It is also not immune to weak information, parameter boundaries, model misspecification, or numerical error. In high-dimensional models, the nuisance optimization may be expensive or difficult, and the resulting profile is only as trustworthy as those constrained fits. The Wald confidence interval and its failure modes are relevant context, but a profile calculation should not be described as a confidence interval until an appropriate inferential rule is applied.

Common Mistake

Wrong reasoning: “I fixed the parameter of interest and evaluated the likelihood while leaving every other parameter at its full-fit estimate, so I have the profile likelihood.” This is a likelihood slice, not a profile. Correction: at every fixed value, reoptimize all parameters that remain free under the model constraints. Otherwise, the computed likelihood can be lower than the maximum attainable under that fixed value and misrepresent the profile.

A second practical mistake is treating an optimizer's success flag as proof that the constrained maximum was found. A success flag reports a stopping condition, not a global-optimality guarantee. Compare starts, inspect gradients and bounds, and verify that the fixed parameter was held exactly at the requested value.

Check Your Understanding

Answer each question, then compare with the brief solutions.

For $\theta=(\psi,\lambda_1,\lambda_2)$, what must be optimized to calculate one profile-likelihood value at $\psi=\psi_0$?
In a Gaussian model with known standard deviations, what quantity can be minimized instead of maximizing the log-likelihood?
A computed profile point has a lower log-likelihood than the unrestricted fit, and the fixed grid contains the unrestricted estimate. Name two checks to perform.
Why is fixing nuisance parameters at their full-fit estimates not the profile likelihood?

Solutions. (1) Maximize $L(\psi_0,\lambda_1,\lambda_2)$ over both nuisance parameters, subject to their constraints. (2) Minimize the weighted sum $\sum_i[(y_i-f_i)/\sigma_i]^2$, since the remaining log-likelihood terms are constant in the parameters. (3) Check the constrained optimizer with multiple starts and inspect whether the fixed value, parameter bounds, likelihood implementation, and numerical settings match the full fit. (4) A profile allows nuisance parameters to adjust to maximize the likelihood at each fixed target value; holding them fixed computes only a likelihood slice.