Tutorials › Quantitative Systems Pharmacology (QSP) › Estimating Interaction Parameters from Factorial Dose Matrices

Combination Therapy · Tutorial 512 of 799

Estimating Interaction Parameters from Factorial Dose Matrices

Fit and evaluate interaction parameters from crossed single-agent and combination dose observations, while keeping the response scale, error model, and limits of the estimate explicit.

Advanced 12 min read

What You'll Learn

  • Write a factorial dose-response model with main effects and an interaction term.
  • Estimate interaction as a difference-in-differences in a two-level design.
  • Fit linear-in-parameter models by least squares and interpret the coefficient uncertainty.
  • Use replicate variation and residual patterns to assess whether the observation model is credible.
  • Calculate a weighted interaction contrast when cell means have different measurement precision.
  • Distinguish a design that estimates coefficients from one that also supports uncertainty assessment.

From a Dose Matrix to an Interaction Estimate

A combination experiment is most informative when it measures the two drugs separately as well as together across a grid of dose pairs. The single-agent arms anchor the response expected from each drug on its own; the combination cells show whether observed responses depart from the chosen reference model. Estimation therefore depends not just on the number of combination wells, but on the doses represented, the response scale, and the observation-error model.

For a first, transparent model, let \(y_{ijr}\) be the measured response at dose levels \(u_i\) and \(v_j\), with replicate \(r\). The doses \(u\) and \(v\) are expressed in stated, fixed units or are normalized by prespecified reference doses. Consider the model

$$ y_{ijr}=\beta_0+\beta_A u_i+\beta_B v_j+\gamma u_i v_j+\epsilon_{ijr}. $$

Here \(\beta_0\) is the expected response at zero dose of both drugs, \(\beta_A\) and \(\beta_B\) are the linear single-agent slopes on this response scale, and \(\gamma\) is the coefficient of the product term. The interaction contribution at a particular dose pair is \(\gamma u_i v_j\); \(\gamma\) itself is the interaction change per dose-unit product. Thus its units are response units divided by the units of \(u\) and \(v\). A positive \(\gamma\) raises the modeled response relative to the additive linear model when both doses are positive; whether that is called synergy depends on what response means and which reference null is appropriate.

Model and observation assumptions: The conditional mean response follows the stated equation over the dose range studied. For ordinary least squares with its usual standard errors and \(t\)-based intervals, assume errors have mean zero, are independent across experimental units, have a common variance \(\sigma^2\), and are normally distributed. Replicates must be independent experimental observations, not merely repeated readings of the same sample treated as independent. The interaction estimate is conditional on the chosen response scale and model form.

The linear dose terms are a useful local model, not a universal pharmacology law. If single-agent responses are curved, replace \(u\) and \(v\) with prespecified response functions, such as \(f_A(u)\) and \(f_B(v)\), and use an interaction basis such as \(f_A(u)f_B(v)\). The key estimation ideas remain: state the mean model, map every matrix cell to its model predictors, and match the error model to how the observations were obtained.

Estimation and the Factorial Contrast

Collect the observations into a vector \(y\) and form a design matrix \(X\) with one row per observation and columns \(1,u,v,uv\). If the model is linear in its coefficients, least squares chooses the coefficient vector that minimizes the residual sum of squares. When \(X\) has full column rank, the solution is

$$ \hat{\boldsymbol{\beta}}=(X^\mathsf{T}X)^{-1}X^\mathsf{T}y, \qquad \hat{\sigma}^2=\frac{\sum_{k=1}^{n}(y_k-\hat{y}_k)^2}{n-p}, $$

where \(\boldsymbol{\beta}=(\beta_0,\beta_A,\beta_B,\gamma)^\mathsf{T}\), \(p=4\) for this model, and \(n\) is the number of observations. Under the normal, independent, constant-variance error assumptions, the estimated coefficient covariance matrix is \(\hat{\sigma}^2(X^\mathsf{T}X)^{-1}\). Its diagonal entry for \(\gamma\) gives the estimated variance of \(\hat{\gamma}\); the square root is \(SE(\hat{\gamma})\).

A balanced two-level factorial design gives a particularly clear interpretation. At doses \(u=0\) or \(u=1\) and \(v=0\) or \(v=1\), let the four fitted cell means be \(\mu_{00},\mu_{10},\mu_{01},\mu_{11}\). Substitution into the model gives \(\mu_{00}=\beta_0\), \(\mu_{10}=\beta_0+\beta_A\), \(\mu_{01}=\beta_0+\beta_B\), and \(\mu_{11}=\beta_0+\beta_A+\beta_B+\gamma\). Subtracting the single-agent changes from the combination change yields

$$ \gamma=\mu_{11}-\mu_{10}-\mu_{01}+\mu_{00}. $$

This is a difference-in-differences: it removes the baseline and both main effects, leaving the departure from additivity in this model. Replicate observations estimate the cell means and, under the stated assumptions, provide information about their uncertainty and the residual variance.

Worked Example: Estimating an Interaction by Difference-in-Differences

Suppose a hypothetical assay records response in assay units at the four dose combinations \(u,v\in\{0,1\}\). The replicate measurements are: \((0,0)\): 10 and 12; \((1,0)\): 18 and 20; \((0,1)\): 16 and 18; \((1,1)\): 30 and 32. Assume independent, normally distributed errors with common variance and the linear factorial mean model above.

The cell means are \(\bar{y}_{00}=11\), \(\bar{y}_{10}=19\), \(\bar{y}_{01}=17\), and \(\bar{y}_{11}=31\) assay units. The interaction estimate is

$$ \hat{\gamma}=31-19-17+11=6 $$

assay units per dose-unit product. The fitted single-agent coefficients are \(\hat{\beta}_A=19-11=8\) assay units per \(u\)-unit and \(\hat{\beta}_B=17-11=6\) assay units per \(v\)-unit; \(\hat{\beta}_0=11\) assay units.

The four within-cell residual pairs are \((-1,1)\), \((-1,1)\), \((-1,1)\), and \((-1,1)\), so the residual sum of squares is \(SSE=8\). There are \(n=8\) observations and \(p=4\) coefficients, giving \(df=8-4=4\) and \(\hat{\sigma}^2=8/4=2\) assay-units squared. Each cell mean has estimated variance \(2/2=1\), because it averages two independent observations. The contrast has four coefficients, each with squared value 1, so \(SE(\hat{\gamma})=\sqrt{1+1+1+1}=2\) assay units per dose-unit product.

The test statistic for \(H_0:\gamma=0\) is \(t=6/2=3.00\), with \(df=4\). The two-sided \(p\)-value is approximately 0.0399. The 95% \(t\) interval is \(6\pm 2.776(2)\), or approximately [0.448, 11.552] assay units per dose-unit product. Under the model, the estimated interaction is positive; the interval is wide, reflecting the small residual degrees of freedom. This alone does not establish a mechanistic interaction or justify calling the result synergy on other response scales.

Design Support, Residuals, and Uncertainty

A single combination dose pair cannot separate the baseline, two main effects, and product contribution if all are unknown: one mean response provides only one equation for four coefficients. However, it is too strong to conclude that every combination-only design is unable to estimate these coefficients. A set of varied dose pairs can provide enough independent predictor combinations, especially if the baseline is known or fixed from appropriate controls. Single-agent arms are valuable because they directly measure single-agent responses and can validate the reference model, not because every combination-only design is mathematically incapable of fitting coefficients.

Worked Example: A Combination-Only Design With a Known Baseline

Consider the same linear form, but suppose the baseline is known to be \(\beta_0=10\) assay units from a reliable control measurement. Three hypothetical combination observations are \(y=6,4,5\) assay units at \((u,v)=(1,1),(1,2),(2,1)\). Subtracting the known baseline gives the equations

$$ \begin{aligned} \beta_A+\beta_B+\gamma&=-4,\\ \beta_A+2\beta_B+2\gamma&=-6,\\ 2\beta_A+\beta_B+2\gamma&=-5. \end{aligned} $$

The design matrix for these three coefficients has rows \((1,1,1),(1,2,2),(2,1,2)\). Its determinant is \(1(4-2)-1(2-4)+1(1-4)=1\), so the rows are linearly independent and the coefficients have a unique solution. Subtracting the first equation from the second gives \(\beta_B+\gamma=-2\); subtracting the first from the third gives \(\beta_A+\gamma=-1\). Substitution into the first equation gives \(\gamma=1\), then \(\beta_A=-2\) and \(\beta_B=-3\).

Thus, these varied combination dose pairs estimate the three coefficients when the baseline is fixed. The three observations exactly determine three coefficients, however, leaving no residual degrees of freedom to estimate an error variance from this fit. Replication or additional informative observations are needed for an empirical uncertainty assessment. If the baseline were also unknown, three observations would not estimate all four coefficients.

A numerical estimate should be accompanied by checks of both mean structure and error structure. Plot residuals \(e_k=y_k-\hat{y}_k\) against fitted values and against each dose; plot them by dose cell, plate position, and experimental run when those factors could matter. A fan-shaped spread suggests nonconstant variance. Systematic curvature or dose-dependent residual runs suggest that linear single-agent terms are inadequate. Clustering by plate or run challenges independence, while an unusual observation may indicate a data or assay issue rather than a genuine interaction.

Examine uncertainty as well as the point estimate. For the linear model under its assumptions, a confidence interval for \(\gamma\) is \(\hat{\gamma}\pm t_{1-\alpha/2,df}SE(\hat{\gamma})\). It describes values compatible with the data and model under repeated sampling; it is not a probability statement about the fixed parameter. Report the coefficient with units, interval, number of independent experimental units, and the response scale. A narrow interval can still be misleading if the model is wrong or experimental replicates were treated as independent when they were not.

Worked Example: Unequal Precision Across Factorial Cells

Suppose hypothetical replicate data produce cell means \(\bar{y}_{00}=10.2\), \(\bar{y}_{10}=18.0\), \(\bar{y}_{01}=16.5\), and \(\bar{y}_{11}=29.7\) assay units. An independently justified measurement-error model gives standard errors for these means of 0.5, 0.6, 0.4, and 0.8 assay units, respectively. Assume the cell means are independent and these standard errors are known.

The factorial interaction contrast is \(\hat{\gamma}=29.7-18.0-16.5+10.2=5.4\) assay units per dose-unit product. Since the four means are independent, the contrast variance is the sum of their variances weighted by the squared contrast coefficients:

$$ SE(\hat{\gamma})= \sqrt{0.5^2+0.6^2+0.4^2+0.8^2} =\sqrt{1.41} =1.187. $$

Using the normal critical value 1.96 for known standard errors, the 95% interval is \(5.4\pm1.96(1.187)\), or approximately [3.073, 7.727] assay units per dose-unit product. The unequal precision is accounted for in the uncertainty calculation; the point contrast remains the same because it is defined by the four cell means. If the standard errors were estimated from small replicate samples rather than known, a \(t\)-based or model-based uncertainty calculation would be needed instead. Weighting cannot repair a misspecified mean model or dependence among the cell means.

Alternative Fits and Practical Limits

Ordinary least squares is appropriate for independent observations with approximately constant variance and a mean that is linear in the fitted coefficients. Weighted least squares is useful when observation variances are known or credibly estimated and differ across observations; it minimizes a variance-weighted residual sum of squares. Nonlinear least squares is appropriate when the response functions are nonlinear in parameters, such as curved single-agent dose-response functions, but its uncertainty approximation can be less reliable in small or poorly supported designs. A hierarchical model can represent plate-to-plate or experiment-to-experiment variation rather than folding it into independent residual noise.

These methods answer different modeling needs; they do not make the interaction parameter independent of the chosen reference and response scale. A checkerboard summary score may be convenient for screening many dose pairs, but it can obscure the fitted single-agent components, the error structure, and parameter uncertainty. Conversely, a parametric interaction coefficient is interpretable only over the dose range and functional form it represents. Extrapolation beyond the matrix is particularly risky.

Common warning signs include highly correlated fitted coefficients, very large standard errors, estimates that change substantially when one dose level is removed, residual curvature, and different apparent interaction directions across response transformations. Such patterns warrant checking dose coverage, baseline handling, replicate structure, and the mean and error models before making a pharmacological claim. A good estimate is not just a coefficient: it is a coefficient tied to a defined response, model, design, and uncertainty calculation.

Common pitfall: “The combination response is lower than either single-agent response, so the fitted interaction must be synergistic.” This reasoning compares outcomes without accounting for the responses expected under the chosen null model. Estimate the single-agent components and the interaction on a specified scale, then interpret the coefficient relative to that model. A coefficient’s sign is not a scale-free biological label.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

  1. In the linear factorial model, what does \(\gamma\) measure, and what are its units?
  2. Cell means at \((0,0),(1,0),(0,1),(1,1)\) are 5, 9, 8, and 15 assay units. Calculate the interaction contrast.
  3. Why can three varied combination dose pairs estimate three coefficients when the baseline is known, while one combination dose pair cannot?
  4. A residual plot shows increasing residual spread as fitted response increases. Name one assumption in question and one reasonable modeling response.
  5. Why is a positive interaction coefficient not, by itself, proof of mechanistic synergy?

Solutions. (1) \(\gamma\) is the product-term contribution per unit of \(u v\), conditional on the specified mean model; its units are response units divided by the product of dose units. (2) \(15-9-8+5=3\) assay units per dose-unit product. (3) With a fixed baseline, three independent rows of dose predictors can determine the three remaining coefficients; a single cell supplies only one equation. Replication or extra observations are still needed to estimate residual variability reliably. (4) Constant error variance may be violated; consider a justified variance model or transformation and then reassess residuals. (5) The sign is conditional on the response scale and reference model, and does not by itself identify a mechanism or establish a scale-independent synergy claim.

Clinical Biostats