Tutorials › Quantitative Systems Pharmacology (QSP) › Profiling Derived Quantities and Parameter Combinations

Practical Identifiability and Estimability · Tutorial 559 of 799

Profiling Derived Quantities and Parameter Combinations

Learn to construct and interpret likelihood profiles for derived quantities while allowing model parameters to adjust subject to the quantity of interest.

Advanced 11 min read

What You'll Learn

  • Define a profile likelihood for a function of model parameters
  • Distinguish profiling a derived quantity from substituting fitted parameters
  • Construct profiles for a parameter ratio, pathway flux, and pathway fraction
  • Calculate likelihood-ratio intervals and interpret them in context
  • Recognize when boundaries, multiple optima, or weak identification complicate a profile
  • Compare profiling with parameter-wise Wald intervals and delta-method approximations

Why Profile a Derived Quantity?

A QSP model is often used to interpret a quantity that is not itself a fitted parameter. Examples include a ratio of two pathway rates, a flux calculated from a rate and a concentration, or the fraction of total pathway activity carried by one route. Such a quantity can matter scientifically even when its component parameters are strongly correlated or individually imprecise.

A direct calculation from the maximum-likelihood estimates gives a point estimate, but does not by itself describe uncertainty in that quantity. Profiling answers a more targeted question: for each proposed value of the derived quantity, how well can the model fit the data after all parameter values compatible with that proposal have been optimized?

Definition: Let \(\theta\) be a vector of model parameters in a feasible parameter space \(\Theta\), let \(L(\theta;y)\) be the likelihood for data \(y\), and let \(g(\theta)\) be a scalar quantity of interest. The profile likelihood for \(q=g(\theta)\) is $$ L_p(q)=\sup_{\theta\in\Theta:\,g(\theta)=q}L(\theta;y). $$ The profile log-likelihood is \(\ell_p(q)=\log L_p(q)\). The unconstrained maximum log-likelihood is \(\hat{\ell}=\sup_{\theta\in\Theta}\log L(\theta;y)\), and the likelihood-ratio statistic at \(q\) is \(D(q)=2\{\hat{\ell}-\ell_p(q)\}\).

The constraint \(g(\theta)=q\) defines all parameter combinations compatible with the proposed value. The optimization searches over that set, rather than holding the other parameters at their unconstrained estimates. The profile is therefore a measure of the best fit available at each \(q\), not the likelihood obtained by plugging a single fixed parameter vector into \(g\).

Constructing and Interpreting a Profile

In practice, first fit the unconstrained model to obtain \(\hat{\theta}\) and \(\hat{\ell}\). Then choose values of \(q\) across its feasible range. For each value, impose \(g(\theta)=q\), optimize the likelihood over the remaining degrees of freedom, and record \(\ell_p(q)\). The profile reaches its maximum at \(q=g(\hat{\theta})\), provided the unconstrained maximum is attained and the optimization is successful.

A likelihood-ratio confidence set contains the values whose loss of fit relative to the unconstrained maximum is no greater than a chosen cutoff \(c\):

$$ \mathcal{C}=\{q: D(q)\le c\}. $$

For a regular, interior, scalar quantity with adequate information, a common large-sample 95% cutoff is \(c=3.84\), the 95th percentile of a chi-squared distribution with one degree of freedom. The cutoff may instead be calibrated for the model and data setting, especially when small samples, boundaries, or nonregular behavior matter. The interval or set describes values not rejected by this likelihood-ratio procedure at the specified level; it is not a probability distribution for the fixed true value.

A profile can be asymmetric, broad, bounded by the model, or even disconnected. These features may be important: a nearly flat profile indicates little information about the quantity, while a narrow profile indicates that substantially different values impose a detectable loss of fit. The profile concerns uncertainty in \(q\), which need not match the uncertainty in any one component parameter.

Worked Examples

Worked Example: A Ratio of Two Pathway Rates

Suppose a hypothetical model has two positive rate parameters \(k_1\) and \(k_2\), both measured in \(\mathrm{h}^{-1}\). The quantity of interest is their dimensionless ratio \(r=k_1/k_2\). For illustration, suppose the fitted log-likelihood can be expressed using \(x=\log(r)\) and \(z=\log\{k_1k_2/(1\ \mathrm{h}^{-2})\}\) as $$ \ell(x,z)=-\frac{1}{2}\left[\frac{(x-\log 2)^2}{0.25}+\frac{(z-\log 10)^2}{9}\right]. $$

The unconstrained maximum occurs at \(\hat{x}=\log 2\), \(\hat{z}=\log 10\), so the fitted ratio is \(\hat{r}=2\). To profile a proposed ratio \(r\), fix \(x=\log r\) and maximize over \(z\). The second squared term is largest in log-likelihood when \(z=\log 10\), so $$ \ell_p(r)=-\frac{1}{2}\frac{(\log r-\log 2)^2}{0.25}, \qquad D(r)=\frac{(\log r-\log 2)^2}{0.25}. $$

For the proposed ratio \(r=4\), $$ D(4)=\frac{(\log 4-\log 2)^2}{0.25} =\frac{(\log 2)^2}{0.25} =\frac{0.4805}{0.25}=1.922. $$ This is below the illustrative 95% cutoff of 3.84, so \(r=4\) is included in the corresponding likelihood-ratio confidence set.

The set is found by solving \(D(r)\le3.84\): $$ |\log r-\log 2|\le\sqrt{3.84(0.25)}=0.9798. $$ Exponentiating both limits gives \(r\in[2e^{-0.9798},2e^{0.9798}]\), or approximately \([0.751,5.33]\). This wide interval is on the ratio scale and is asymmetric around 2 in ordinary arithmetic, even though it is symmetric on the log-ratio scale.

Worked Example: Profiling a Pathway Flux

Consider a hypothetical pathway flux \(J=kC\), where \(k\) is in \(\mathrm{h}^{-1}\) and \(C\) is in \(\mathrm{nM}\), making \(J\) measured in \(\mathrm{nM\,h}^{-1}\). Define dimensionless parameters \(A=k/(1\ \mathrm{h}^{-1})\) and \(B=C/(1\ \mathrm{nM})\); then the numerical value of \(J\) in \(\mathrm{nM\,h}^{-1}\) is \(q=AB\). Write \(a=\log A\) and \(b=\log B\). Suppose the hypothetical log-likelihood is $$ \ell(a,b)=-\frac{1}{2}\left[\frac{(a+b-\log 12)^2}{0.04}+\frac{(a-b)^2}{4}\right]. $$

The unconstrained maximum has \(a+b=\log 12\) and \(a-b=0\), giving \(A=B=\sqrt{12}\) and \(\hat{J}=12\ \mathrm{nM\,h}^{-1}\). To profile a proposed flux \(q\), impose \(a+b=\log q\). At this fixed sum, the first squared term is fixed, while the second is maximized in log-likelihood at \(a-b=0\). Thus the compatible parameter values giving the best fit have \(A=B=\sqrt q\), and $$ D(q)=\frac{(\log q-\log 12)^2}{0.04}. $$

Using \(c=3.84\), the interval satisfies $$ |\log q-\log 12|\le\sqrt{3.84(0.04)}=0.3919. $$ Therefore \(q\) ranges from \(12e^{-0.3919}=8.11\) to \(12e^{0.3919}=17.76\), giving an illustrative interval of \([8.11,17.8]\ \mathrm{nM\,h}^{-1}\).

The example also shows why the derived quantity can be more informative than its components. The likelihood penalizes changes in \(a-b\) relatively little, so many rate-and-concentration pairs can fit similarly, but it penalizes changes in \(a+b=\log(AB)\) more strongly. The flux can consequently have a more focused profile than either component considered alone.

Worked Example: A Fraction of Total Pathway Activity

Suppose two positive pathway contributions have a fraction of interest \(q=J_1/(J_1+J_2)\), with \(0<q<1\). Let \(p=\log\{q/(1-q)\}\) be its logit and let \(s\) denote the log of a suitably scaled total contribution. Assume the hypothetical log-likelihood is $$ \ell(p,s)=-\frac{1}{2}\left[\frac{p^2}{0.09}+(s-\log 10)^2\right]. $$

The unconstrained estimates are \(\hat{p}=0\) and \(\hat{s}=\log 10\), so \(\hat{q}=0.5\). For a fixed fraction \(q\), \(p=\log\{q/(1-q)\}\). Maximizing over \(s\) gives \(s=\log 10\), and hence $$ D(q)=\frac{\left[\log\{q/(1-q)\}\right]^2}{0.09}. $$

For \(q=0.60\), the logit is \(\log(0.60/0.40)=\log(1.5)=0.4055\), so $$ D(0.60)=\frac{0.4055^2}{0.09}=\frac{0.1644}{0.09}=1.827. $$ This value is below 3.84. The illustrative 95% profile set is obtained from \(|p|\le\sqrt{3.84(0.09)}=0.5879\). Applying the logistic transformation to the two endpoints gives \(q\in[0.357,0.643]\). These endpoints are on the fraction scale; they are not obtained by adding and subtracting a constant from \(\hat{q}\).

Why the Profile Calculation Works

For a fixed proposed \(q\), the constrained optimizer chooses the nuisance coordinates to make the fit as good as that constraint permits. In the flux example, define \(u=a+b\) and \(v=a-b\). The flux constraint fixes \(u=\log q\), leaving \(v\) free. Substitution into the log-likelihood gives $$ \ell=-\frac{1}{2}\left[\frac{(\log q-\log 12)^2}{0.04}+\frac{v^2}{4}\right]. $$

The first term cannot change while \(q\) is fixed. Since \(v^2\ge0\), the largest possible log-likelihood occurs at \(v=0\). Subtracting this constrained maximum from the unconstrained maximum, which is zero in the displayed parameterization, and multiplying by 2 yields $$ D(q)=\frac{(\log q-\log 12)^2}{0.04}. $$ This argument is the core of profiling: optimize along the constraint set, then compare the best constrained fit with the best fit without the constraint.

Diagnostics, Alternatives, and Limitations

A valid profile requires a well-defined likelihood, a correctly specified feasible parameter space, and successful constrained optimization at each grid value. The constraint must represent the intended scientific quantity, including units and any transformations. For example, a ratio of two rates is dimensionless, whereas a rate multiplied by a concentration has flux units.

In numerical work, inspect the fitted parameters at each profile point as well as the profile curve. A sudden jump in nuisance estimates, abrupt discontinuity, or a profile that depends strongly on the starting values can signal competing local optima or failed optimization. Repeating fits from multiple starting values and checking that \(g(\theta)=q\) is satisfied helps distinguish a genuine shape from a numerical artifact. Also check the profile at its maximum: it should agree with the unconstrained fit up to numerical tolerance.

Profiling has useful strengths and important limits. It accommodates nonlinear transformations and parameter compensation without relying on a symmetric approximation around \(\hat{\theta}\). However, a profile may be computationally expensive, and its confidence set can be unreliable if the likelihood-ratio cutoff is inappropriate. At boundaries, with weakly identified or non-smooth quantities, or when the profile has multiple separated regions, standard regularity conditions may fail; use a calibration appropriate to that setting rather than treating the usual cutoff as automatic.

ApproachWhat it does for a derived quantityMain caution
Profile likelihoodFixes the quantity and re-optimizes all compatible parametersRequires repeated constrained optimization and suitable calibration
Plug-in estimateEvaluates \(g(\hat{\theta})\) at the fitted parameter vectorProvides a point estimate, not an uncertainty assessment
Delta-method approximationApproximates uncertainty using a local linearization near \(\hat{\theta}\)Can be misleading for nonlinear transformations or non-symmetric uncertainty
Reparameterized fitUses the quantity as an explicit model coordinate, if feasibleStill requires accounting for the other parameters and preserving the model constraints

A common mistake is to compute separate intervals for two parameters and combine their endpoints to make an interval for a ratio or product. That does not generally produce a confidence interval for the derived quantity: it ignores the joint likelihood and parameter dependence, and the resulting bounds need not have the intended coverage. The correction is to profile the quantity directly, or use a justified approximation that accounts for the joint uncertainty.

Key takeaway: A profile for \(q=g(\theta)\) is built by fixing \(q\), optimizing over every compatible parameter combination, and comparing that best constrained fit with the unconstrained maximum. It directly targets uncertainty in the scientific quantity of interest.

Check Your Understanding

Use the definitions and calculations above to answer the questions. Solutions follow each question.

  1. If the unconstrained maximum log-likelihood is \(-20.0\) and the constrained maximum at \(q_0\) is \(-22.1\), calculate \(D(q_0)\). With cutoff 3.84, is \(q_0\) included?

    Solution: \(D(q_0)=2\{-20.0-(-22.1)\}=4.20\). Since \(4.20>3.84\), \(q_0\) is not included in that likelihood-ratio set.

  2. For the flux profile in the tutorial, calculate \(D(6)\) using \(D(q)=(\log q-\log12)^2/0.04\). Is \(6\ \mathrm{nM\,h}^{-1}\) included under cutoff 3.84?

    Solution: \(\log6-\log12=-\log2=-0.6931\), so \(D(6)=0.6931^2/0.04=12.01\). Since \(12.01>3.84\), the proposed flux is not included.

  3. Why must the nuisance parameters be re-optimized at each fixed value of \(q\)?

    Solution: The profile asks for the best fit available subject to \(g(\theta)=q\). Holding nuisance parameters at their unconstrained estimates may produce a poorer fit than another compatible parameter combination, so it does not in general calculate \(L_p(q)\).

  4. For a fraction with logit profile \(D(q)=p^2/0.09\), where \(p=\log\{q/(1-q)\}\), calculate \(D(0.40)\). Is it included with cutoff 3.84?

    Solution: The logit is \(\log(0.40/0.60)=-0.4055\). Thus \(D(0.40)=0.4055^2/0.09=1.827\), which is below 3.84, so the fraction is included.

Clinical Biostats