Tutorials › Pharmacokinetic-Pharmacodynamic (PK/PD) Modeling › Profile Likelihood for Indirect Model Parameters

Indirect Response Models · Tutorial 140 of 598

Profile Likelihood for Indirect Model Parameters

You will learn how to profile indirect-model parameters, interpret likelihood-ratio intervals, and recognize when data leave a parameter direction unidentifiable.

Advanced 12 min read

What You'll Learn

  • Define a profile likelihood by maximizing over nuisance parameters at each fixed parameter value.
  • Calculate a likelihood-ratio statistic and a profile-based interval for kout.
  • Interpret a broad but bounded Emax profile as finite-sample uncertainty.
  • Recognize an exactly flat IC50–Emax profile direction as structural non-identifiability.
  • Check profile results for boundaries, optimization failures, and sensitivity to assumptions.
  • Compare profile intervals with Wald intervals and joint likelihood contours.

Why Profile an Indirect-Model Parameter?

A fitted indirect response model can reproduce the observed response well while leaving one or more parameter values uncertain. This is especially relevant for \(k_{out}\), \(E_{max}\), and \(IC_{50}\): response timing and effect magnitude may provide different amounts of information about each, and the parameters can compensate for one another. A standard error from a local approximation may not reveal a long, asymmetric, or nearly flat direction in the likelihood.

A profile likelihood makes that direction visible. It asks: if one parameter is fixed at a particular value, how well can the remaining parameters fit the data? Repeating that calculation across plausible fixed values shows both the best-supported region and whether the data genuinely constrain the parameter.

Definition: Let \(\psi\) be the parameter of interest and \(\lambda\) denote all other model parameters. The profile log-likelihood is \( \ell_p(\psi)=\max_{\lambda}\ell(\psi,\lambda) \). At each fixed value of \(\psi\), the nuisance parameters \(\lambda\) are re-estimated to maximize the log-likelihood. The profile likelihood is the likelihood corresponding to this maximized log-likelihood.

For an inhibitory indirect response model, one possible specification is \( \frac{dR}{dt}=k_{in}\left[1-\frac{E_{max}C(t)}{IC_{50}+C(t)}\right]-k_{out}R(t) \). Here \(R(t)\) is the response, \(C(t)\) is the drug concentration, and \(k_{in}\) and \(k_{out}\) are input and loss rates. In this example, \(E_{max}\) is a dimensionless maximum fractional inhibition and \(IC_{50}\) has the same concentration units as \(C(t)\). The parameter bounds must keep the input rate meaningful; for this form, one commonly imposes \(0\leq E_{max}\leq1\), \(IC_{50}>0\), and \(k_{out}>0\). Other indirect models can use a different effect function or parameterization, so the likelihood must be profiled under the specific model being fitted.

Suppose observations \(Y_i\) at times \(t_i\) are conditionally independent and normally distributed around model predictions \(R(t_i;\theta)\), with a common known residual standard deviation \(\sigma\). Then \(Y_i\sim N(R(t_i;\theta),\sigma^2)\), and the log-likelihood, up to a constant that does not depend on the parameters, is \( \ell(\theta)=-\frac{1}{2\sigma^2}\sum_i[Y_i-R(t_i;\theta)]^2 \). Thus maximizing the likelihood is equivalent to minimizing the sum of squared residuals when \(\sigma\) is fixed. If \(\sigma\) is unknown, it must also be estimated, including in each constrained fit.

Likelihood-Ratio Cutoffs and Their Justification

Let \(\hat{\theta}\) denote the unrestricted maximum-likelihood estimate and \(\hat{\psi}\) its component of interest. The profile likelihood-ratio statistic at a candidate value \(\psi\) is \( \Lambda(\psi)=2[\ell(\hat{\theta})-\ell_p(\psi)] \). Small values mean that fixing \(\psi\) causes little loss of fit; large values mean the constrained fit is substantially worse than the unrestricted fit.

For one parameter of interest, a common approximate 95% likelihood-ratio confidence set contains values satisfying \(\Lambda(\psi)\leq3.84\), the 95th percentile of a chi-square distribution with one degree of freedom. The justification is Wilks' likelihood-ratio result: under regularity conditions, with an identifiable interior parameter and adequate information, twice the log-likelihood difference between the unrestricted and one-constraint fits is approximately chi-square with one degree of freedom. The cutoff is an approximation, not a guarantee. Boundary estimates, small samples, weak identifiability, or a flat likelihood can make the approximation unreliable.

For fixed \(\sigma\), the calculation is particularly transparent. Writing \(SSE_{min}\) for the unrestricted residual sum of squares and \(SSE_p(\psi)\) for the constrained, minimized sum at fixed \(\psi\), \( \Lambda(\psi)=\frac{SSE_p(\psi)-SSE_{min}}{\sigma^2} \). The denominator is \(\sigma^2\), not \(\sigma\), because the log-likelihood is proportional to squared residuals divided by \(2\sigma^2\).

Worked Example: Profiling the Turnover Rate

Consider hypothetical response measurements from 12 independent sampling times. A Gaussian residual model with known \(\sigma=0.05\) response units is used. The unrestricted fit gives \(SSE_{min}=0.030\) response-units squared and \(\hat{k}_{out}=0.16\ \mathrm{h}^{-1}\). At each fixed \(k_{out}\), the other model parameters are re-optimized. The resulting profile values are:

Fixed \(k_{out}\) (\(\mathrm{h}^{-1}\))Profile SSE (response-units squared)\(\Lambda=(SSE_p-SSE_{min})/0.0025\)
0.120.0404.00
0.140.0352.00
0.160.0300.00
0.180.0341.60
0.200.0404.00

For \(k_{out}=0.12\ \mathrm{h}^{-1}\), for example, \( \Lambda=(0.040-0.030)/(0.05^2)=0.010/0.0025=4.00 \). This exceeds 3.84, so that value lies just outside the approximate 95% profile interval. At \(0.14\ \mathrm{h}^{-1}\), the statistic is \(0.005/0.0025=2.00\), which is below the cutoff.

Linear interpolation of the profile statistic between the displayed points gives a lower crossing of about \(0.122\ \mathrm{h}^{-1}\): \(0.12+(4.00-3.84)(0.02)/(4.00-2.00)=0.1216\ \mathrm{h}^{-1}\). The upper crossing is about \(0.199\ \mathrm{h}^{-1}\): \(0.18+(3.84-1.60)(0.02)/(4.00-1.60)=0.1987\ \mathrm{h}^{-1}\). The approximate interval is therefore \(0.122\) to \(0.199\ \mathrm{h}^{-1}\). This is an interval for the turnover rate under the stated model and error assumptions, not a probability statement about the fixed parameter after seeing the data.

Reading Profiles for Effect Parameters

To profile \(E_{max}\), fix it at each candidate value and re-optimize \(k_{out}\), \(IC_{50}\), and any other estimated parameters. The same procedure applies to \(IC_{50}\): hold it fixed and re-optimize the remaining parameters. Re-optimization is essential. Comparing fits while leaving nuisance parameters at their unrestricted estimates is not a profile likelihood and can exaggerate how much information the data contain about the fixed parameter.

Worked Example: A Broad but Finite-Sample Emax Profile

Suppose a hypothetical study has response observations at several concentrations, but the data provide only limited separation of the maximum effect from potency. Assume again that \(n=12\), \(\sigma=0.05\), and \(SSE_{min}=0.030\). The following profile SSE values result after re-optimizing \(IC_{50}\), \(k_{out}\), and the other nuisance parameters at each fixed \(E_{max}\):

Fixed \(E_{max}\)Profile SSE\(\Lambda=(SSE_p-0.030)/0.0025\)
0.200.0424.80
0.400.0320.80
0.600.03050.20
0.800.0300.00
0.900.0310.40
1.000.0404.00

At \(E_{max}=0.60\), \(\Lambda=(0.0305-0.030)/0.0025=0.20\), so this value fits almost as well as the maximum-likelihood estimate \(0.80\). At \(E_{max}=0.20\), the statistic is \(0.012/0.0025=4.80\), which is above 3.84. Interpolating the lower crossing between 0.20 and 0.40 gives \(0.20+(4.80-3.84)(0.20)/(4.80-0.80)=0.248\). Interpolating the upper crossing between 0.90 and 1.00 gives \(0.90+(3.84-0.40)(0.10)/(4.00-0.40)=0.996\). The approximate profile interval is about \(0.248\) to \(0.996\).

This broad, yet crossing, profile indicates substantial finite-sample uncertainty: a wide range of \(E_{max}\) values remains compatible with the data under the likelihood-ratio criterion. It is not an exactly flat direction. In a real analysis, the likelihood should be evaluated at a sufficiently fine grid, and the reported crossings should be checked with direct constrained fits rather than relying on coarse linear interpolation.

Worked Example: An Exactly Flat Emax–IC50 Direction

Now consider a hypothetical steady-state observation at only one concentration, \(C=2\ \mathrm{mg/L}\). Suppose the observed fractional inhibition is \(q=0.30\). For the inhibitory Emax function, \( q=\frac{E_{max}C}{IC_{50}+C} \). Solving this equation for \(IC_{50}\) gives \( IC_{50}=C\left(\frac{E_{max}}{q}-1\right) \).

For \(E_{max}=0.50\), the matching potency is \(IC_{50}=2(0.50/0.30-1)=1.333\ \mathrm{mg/L}\). For \(E_{max}=0.80\), it is \(IC_{50}=2(0.80/0.30-1)=3.333\ \mathrm{mg/L}\). For \(E_{max}=1.00\), it is \(IC_{50}=2(1.00/0.30-1)=4.667\ \mathrm{mg/L}\). Each pair predicts the same inhibition, \(q=0.30\), at the observed concentration. If the likelihood depends on the predicted response at this sole concentration, all these pairs give the same likelihood.

Consequently, when profiling \(E_{max}\) and re-optimizing \(IC_{50}\), the profile is exactly flat over the feasible range of matching pairs (here \(0.30<E_{max}\leq1\), with the lower endpoint approached as \(IC_{50}\) approaches zero). Increasing measurement precision at that same concentration would not separate the two parameters: the problem is structural for this design, not simply a wide finite-sample interval. Additional information that changes the response predictions differently across candidate pairs would be needed to distinguish them.

Diagnostics, Alternatives, and Limitations

A profile is informative only if its constrained optimizations are reliable. Check that each fit converges, that nuisance parameters are re-optimized from multiple starting values when the surface is difficult, and that the best constrained likelihood never exceeds the unrestricted maximum. Extend the parameter grid until the profile clearly crosses the chosen cutoff or reaches a meaningful model boundary. A profile that appears flat only because the grid is too narrow is not evidence of non-identifiability.

A finite-sample profile generally has a maximum and eventually rises as a parameter moves away from it; the uncertainty interval is the region below the selected likelihood-ratio cutoff. An exactly flat profile over a continuum indicates that adjusting nuisance parameters can preserve the likelihood exactly, as in the single-concentration example. A profile that is very shallow but eventually crosses the cutoff instead indicates weak information and finite-sample uncertainty. A profile that does not cross before the allowed parameter boundary yields a one-sided or boundary-limited result, not a two-sided interval with an invented endpoint.

A Wald interval uses a local quadratic approximation around the estimate, often of the form estimate plus or minus a standard error times a normal quantile. It can be convenient when the likelihood is approximately symmetric and well-curved, but it may give misleading ranges for positive, bounded, or highly correlated parameters. A profile interval uses the fitted likelihood away from the maximum and can reflect asymmetry and nuisance-parameter compensation. Joint likelihood contours are useful for displaying dependence between two parameters, such as \(E_{max}\) and \(IC_{50}\), whereas a one-parameter profile summarizes support for one parameter after optimizing the rest. Neither approach repairs a model or design that does not identify the parameters.

The cutoff based on the chi-square approximation also has limitations. If a parameter is at a boundary, the usual reference distribution may not apply. With weak information or a small sample, likelihood-ratio coverage can differ from its nominal level. The profile itself is conditional on the chosen indirect model, observation-error distribution, and treatment of residual variability. It describes uncertainty within those assumptions; it does not establish that the mechanism is correct.

Common pitfall: A narrow standard error from the optimizer does not prove that \(k_{out}\), \(E_{max}\), or \(IC_{50}\) is well identified. A local curvature calculation can miss a distant flat ridge or a boundary. Inspect the profile, verify the constrained fits, and distinguish a finite likelihood-ratio crossing from a profile that remains flat or fails to cross.

Check Your Understanding

Answer each question, then compare with the solutions below.

  1. What is re-optimized when profiling \(IC_{50}\), and why?
  2. With known \(\sigma=0.04\), \(SSE_{min}=0.020\), and constrained \(SSE_p=0.0248\), calculate the likelihood-ratio statistic. Does the candidate lie inside the approximate 95% profile set?
  3. A profile is nearly flat over a wide range but crosses the 3.84 cutoff at both ends. Is this exact non-identifiability or finite-sample uncertainty?
  4. At a single concentration, several \(E_{max}\)–\(IC_{50}\) pairs give exactly the same predicted inhibition. Would collecting a more precise response at that same concentration necessarily separate the parameters?

Solutions. (1) Fix \(IC_{50}\) at each candidate value and re-optimize every other estimated parameter, because the profile measures the best likelihood attainable with that value imposed. (2) \(\Lambda=(0.0248-0.020)/(0.04^2)=0.0048/0.0016=3.00\). Since \(3.00<3.84\), the candidate is inside the approximate 95% profile set. (3) It is finite-sample uncertainty: the likelihood eventually excludes values beyond the crossings, though the broad profile indicates weak information. (4) No. If those pairs yield identical predictions at that concentration, greater precision there does not break the exact parameter trade-off.

Clinical Biostats