Likelihood-Ratio Tests for Nested PK/PD Models

Why Test Nested PK/PD Models?

A PK/PD analysis may compare a simpler exposure-response structure with a more flexible one. For example, a model with no linear exposure effect can be compared with a model that adds an exposure slope. If the simpler model is a special case of the more flexible model, the models are nested. A likelihood-ratio test (LRT) quantifies how much the best achievable likelihood improves when the additional parameter or parameters are allowed.

The test addresses a specific question: is the improvement in fit larger than would typically be expected under the simpler model, given the assumptions used to calibrate the test? It does not establish that a statistically significant added effect is clinically important, nor does a nonsignificant result prove that the simpler model is true.

Definition: Let $M_0$ be a reduced model and $M_1$ a full model. They are nested when every probability distribution allowed by $M_0$ is also allowed by $M_1$, usually by fixing one or more parameters of $M_1$ to specified values. Write $\ell_0$ and $\ell_1$ for their maximized log likelihoods, fitted to the same observed data under comparable likelihood formulations. The likelihood-ratio statistic is $$ \Lambda = 2(\ell_1-\ell_0). $$ Because the full model contains the reduced model, $\ell_1\geq \ell_0$ and $\Lambda\geq 0$, apart from numerical fitting error.

For a test of the reduced model, the null hypothesis $H_0$ states that the restrictions defining $M_0$ hold; the alternative $H_a$ allows the additional parameter values in $M_1$. If the models differ by $q$ freely varying parameters, $q$ is the nominal degrees of freedom for the regular LRT. The number of parameters means the difference in the dimensions of the models, not the number of observations or the number of model labels.

Derivation and the Regular Chi-Square Test

The likelihood $L$ is the probability density or probability mass assigned to the observed data as a function of model parameters. Maximizing it separately under the two models gives $L_0$ and $L_1$, with $L_1\geq L_0$. The likelihood ratio in favor of the reduced model relative to the full model is $L_0/L_1$, a number between zero and one. Taking minus twice its logarithm gives

-2\log\left(\frac{L_0}{L_1}\right) = -2(\ell_0-\ell_1) = 2(\ell_1-\ell_0) = \Lambda.

The statistic is larger when the full model improves the maximized likelihood more. Under the regularity conditions of Wilks' theorem, if $H_0$ is true and the sample size is sufficiently large, $\Lambda$ has approximately a chi-square distribution with $q$ degrees of freedom. The p-value is the upper-tail probability $P(\chi^2_q\geq\Lambda_{\mathrm{obs}})$, calculated under that null reference distribution.

The regular approximation requires more than nesting. The true null parameter should be an interior point of the full model's parameter space; parameters should be identifiable; the likelihood should be sufficiently smooth; and the information about the parameters should be nondegenerate. The data-generating likelihood must also be adequately specified for the reference distribution to be meaningful. For PK/PD data, the likelihood and fitting method must account for the relevant observation distribution and dependence structure. A nominal chi-square test is not made valid merely by having two fitted models and a difference in parameter counts.

Worked Example: Adding a Linear Exposure Effect

Hypothetical data. Suppose six independent participants have one PD response measured at a fixed time, with their exposure values coded $x=0,0,1,1,2,2$. The corresponding responses are $y=1,2,2,3,4,5$ response units. Compare a constant-mean model $M_0$, $Y_i\sim N(\mu,\sigma^2)$, with $M_1$, $Y_i\sim N(\beta_0+\beta_1x_i,\sigma^2)$. Both models estimate the same residual variance. Assume independent normal residuals with a common variance, fixed exposure values, and no unmodeled repeated-measure dependence. The slope $\beta_1$ is unrestricted, so its null value zero is an interior point.

Under $M_0$, the fitted mean is $\bar{y}=17/6=2.833$. Since $\sum y_i^2=59$, its residual sum of squares is $$ RSS_0=\sum y_i^2-\frac{(\sum y_i)^2}{6} =59-\frac{17^2}{6} =\frac{65}{6}=10.833. $$ Under $M_1$, $\bar{x}=1$, $\bar{y}=17/6$, and the fitted slope is $$ \hat{\beta}_1= \frac{\sum (x_i-\bar{x})(y_i-\bar{y})} {\sum (x_i-\bar{x})^2} =\frac{6}{4}=1.500. $$ The intercept is $\hat{\beta}_0=\bar{y}-\hat{\beta}_1\bar{x}=1.333$, and the full-model residual sum of squares is $RSS_1=11/6=1.833$. Thus the fitted response increases by 1.500 response units per coded exposure unit; the LRT tests whether the data support a nonzero slope under the model assumptions.

For a normal model with an estimated common variance, the maximized log likelihood is $\ell=-\frac{n}{2}\{\log(2\pi)+1+\log(RSS/n)\}$. The terms common to both models cancel, so $$ \Lambda=n\log\left(\frac{RSS_0}{RSS_1}\right) =6\log\left(\frac{65/6}{11/6}\right) =6\log(65/11) =10.66. $$ The models differ by one freely varying parameter, the slope, so the regular reference distribution is $\chi^2_1$. The upper-tail p-value is approximately $0.0011$. Under the stated assumptions, this is evidence against a zero linear exposure slope. It is not a probability that $H_0$ is true, and it does not by itself establish that a linear exposure-response model is scientifically adequate.

When the Null Is on a Boundary

PK/PD parameters are often constrained by scientific meaning or the model's mathematical form. A drug-effect magnitude may be constrained to be nonnegative. Then testing a zero effect places the null value at the edge, or boundary, of the allowed parameter space rather than in its interior. The regular chi-square result does not automatically apply.

In a common special case, a single parameter $b\geq 0$ is tested at $H_0:b=0$; all nuisance parameters remain identifiable and satisfy regular conditions, and the likelihood has positive information in the relevant direction. The limiting null distribution is then a mixture: with probability one-half the statistic is zero, and with probability one-half it follows $\chi^2_1$. For a positive observed statistic, the asymptotic p-value is one-half the upper-tail $\chi^2_1$ probability. This result is specific to the stated conditions; it is not a universal correction for every constrained PK/PD parameter.

Worked Example: A Nonnegative Drug-Effect Magnitude

Hypothetical data. Consider four independent placebo-adjusted PD responses, in response units, $y=0,1,2,3$. For illustration, suppose their model is $Y_i\sim N(b,4)$, where the residual variance of $4$ response-units squared is known and the drug-effect magnitude is constrained to $b\geq0$. The reduced model fixes $b=0$. Assume the normal model and independence are appropriate, and that there are no additional unidentified parameters.

The sample mean is $\bar{y}=6/4=1.500$, so the constrained full-model estimate is $\hat b=1.500$. With known variance $\sigma^2=4$, the log-likelihood improvement is $$ \ell_1-\ell_0 =\frac{\sum y_i^2-\sum(y_i-\bar{y})^2}{2\sigma^2} =\frac{14-5}{8} =1.125. $$ Therefore, $\Lambda=2(1.125)=2.250$. Because $b=0$ is the boundary of $b\geq0$, use the half-$\chi^2_1$, half-point-mass mixture under the conditions above. The p-value is $$ p=\frac{1}{2}P(\chi^2_1\geq2.250) \approx \frac{1}{2}(0.1336) =0.0668. $$ At a significance level of $0.05$, this test does not provide sufficient evidence against the zero-effect model. The data and model assumptions remain important: the mixture reference is justified here by the simple one-boundary-parameter setup, not simply because the model concerns a drug effect.

Zero Emax and Non-Identifiability

A particularly important warning arises with an Emax exposure-response model. For example, a mean response may contain the term $E_{\max}C/(EC_{50}+C)$, where $C$ is exposure. If the null hypothesis sets $E_{\max}=0$, the mean no longer depends on $EC_{50}$. Under the null, $EC_{50}$ cannot be identified from the response likelihood. This violates a key regularity condition behind both the ordinary Wilks result and the simple one-parameter boundary mixture.

Consequently, neither “one added parameter means $\chi^2_1$” nor “the effect is nonnegative, so halve the $\chi^2_1$ p-value” is generally justified for a zero-Emax test. The null distribution can depend on the exposure design and on how the unidentified parameter is searched or constrained. A parametric bootstrap can calibrate the test for the specified design and fitting procedure, provided the null model is adequate and the simulation and refitting steps faithfully reproduce the analysis.

Worked Example: Calibrating a Zero-Emax Test by Bootstrap

Hypothetical data. Suppose a study has 24 participants assigned across placebo and three exposure levels, with a continuous PD response. The reduced model has no Emax drug effect; the full model adds $E_{\max}\geq0$ and $EC_{50}>0$. Assume the chosen residual likelihood is appropriate, and both models are fitted to the same responses with the same baseline and residual structure. Here, $EC_{50}$ is unidentified when $E_{\max}=0$, so use a parametric-bootstrap reference rather than assuming a standard chi-square mixture.

Suppose the fitted log likelihoods are $\ell_0=-84.2$ and $\ell_1=-81.4$. The observed statistic is $$ \Lambda_{\mathrm{obs}}=2\{\ell_1-\ell_0\} =2\{-81.4-(-84.2)\} =5.600. $$ To calibrate it, simulate 2,000 data sets from the fitted null model at the observed design. For each simulated data set, refit both models using the same parameter constraints, exposure search, and optimization procedure, then calculate its LRT statistic. Suppose 74 simulated statistics are at least 5.600. With the common finite-simulation correction, $$ \hat p=\frac{74+1}{2000+1} =\frac{75}{2001} =0.0375. $$ This estimated bootstrap p-value is evidence against the null model at the 0.05 level, conditional on the null model and bootstrap procedure. The correction avoids reporting a zero p-value when no simulated statistic exceeds the observed one; it does not remove uncertainty from using a finite number of simulations. Increasing the number of simulations improves the numerical precision of the estimated tail probability.

Model Comparability, Diagnostics, and Alternatives

An LRT is meaningful only when the likelihoods are genuinely comparable. Both fits must use the same observations, response scale, and likelihood formulation. For a fixed-effects comparison in a mixed-effects model, likelihoods based on maximum likelihood are generally required; restricted maximum likelihood likelihoods are not directly comparable when the fixed-effects structures differ. The random-effects and residual structures should be handled consistently unless their change is itself the defined model comparison. AIC or other penalized criteria are different tools for model comparison and do not supply the LRT's null p-value.

Before interpreting the statistic, check that the reduced model is truly obtainable by restricting the full model, that both optimizations converged to credible maxima, and that the full-model log likelihood is not lower due to a failed fit. In nonlinear PK/PD models, inspect parameter estimates, standard errors or curvature, boundary estimates, and sensitivity to starting values. Near-flat likelihoods, strong parameter correlations, or estimates that change markedly across starts can signal weak identification. The earlier material on PD diagnostics helps assess whether the candidate model's residual and predictive behavior is credible; an LRT cannot repair a misspecified likelihood or inadequate candidate structure.

Wald tests assess an estimate relative to an estimated standard error, while score tests assess the likelihood's local slope at the null. Under regular conditions, these approaches can be asymptotically related, but they can disagree in finite samples, especially with nonlinear models, weak identification, or boundary constraints. Bootstrap calibration is useful when the reference distribution is nonstandard and the data-generating assumptions can be simulated. For non-nested models, an LRT is not appropriate because neither model is a restricted case of the other.

Common pitfall: “The full model adds one parameter, so compare the LRT with $\chi^2_1$.” Parameter counting gives the regular degrees of freedom only when the null lies in the interior and the parameters are identifiable with the required regularity. For a constrained zero effect, check whether a boundary mixture applies; for zero Emax with unidentified $EC_{50}$, use a suitable nonstandard calibration such as a parametric bootstrap.

LRTs are useful because they give a transparent, likelihood-based test of a precisely stated restriction. Their limitations are equally important: a small p-value does not establish clinical relevance, a large p-value does not show equivalence, and the test compares only the restrictions represented by the chosen nested pair. Report the statistic, its reference distribution and degrees of freedom when applicable, the p-value, the model assumptions, and the scientific meaning and units of the added effect.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

A full model has maximized log likelihood $-46.1$, and its nested reduced model has maximized log likelihood $-48.6$. Calculate the LRT statistic.
A model adds two freely varying, identifiable parameters, and the null values are interior points of their parameter spaces. Under regular conditions, what reference distribution is used?
Why is the half-$\chi^2_1$, half-point-mass reference not automatically valid for a test of zero Emax?
In a bootstrap test, 19 of 999 simulated statistics are at least as large as the observed statistic. Calculate the corrected estimated p-value.
Two models use different observations. Can their maximized likelihoods ordinarily be compared with a standard nested-model LRT? Explain.

Solutions. (1) $\Lambda=2\{-46.1-(-48.6)\}=2(2.5)=5.00$. (2) Use $\chi^2_2$, provided the other regularity and likelihood assumptions also hold. (3) At $E_{\max}=0$, the response no longer depends on $EC_{50}$, so $EC_{50}$ is unidentified under the null; the simple boundary-mixture result assumes the remaining parameters are identifiable and regular. (4) The corrected estimate is $(19+1)/(999+1)=20/1000=0.0200$. (5) No. The fits must use the same observed data and comparable likelihood formulations; a change in observations changes the likelihood comparison and invalidates the ordinary nested-model LRT calibration.