Analysis of multiplex immunoassay data subject to left and right censoring

Elizabeth G. Hill^{1} Hollings Cancer Center Department of Biostatistics, Bioinformatics and Epidemiology Medical University of South Carolina Charleston, SC 29425 (843) 876-1115 fax: (843) 792-4233 hille@musc.edu

1 contact author

1 Introduction

The four- and five-parameter logistic models are examples of nonlinear dose-response curves commonly used for calibration and dose-estimation in multiplex immunoassays (e.g. Bio-Rad’s Bioplex^{\begin{matrix}c \\ ◯ \end{matrix}} system). Let y be the response (e.g. flourescence) at dose x. Then the form of the five-parameter logistic (5PL) model is

y = f(x, β) = d + \frac{a − d}{(1 + (\frac{x}{c})^{b})^{g}} ,  (1)

where β = (a, b, c, d, g)′. The four-parameter logistic model (4PL) is equivalent to the 5PL with g = 1.

The dose-response curve is fit based on replicates of known standards. In my experience, the number of replicates is typically two and the number of standards is eight to ten. The standard concentrations are typically manufacturer-recommended serial dilutions of the analyte of interest. Figure 1 shows the fitted standard curve based on two replicates for each of eight standard concentrations of IL7. The shape of the curve is sigmoidal when the response is graphed as a function of log-dose. The standard curve is used to estimate analyte concentrations for unknowns (e.g. patient urine or serum samples) based on the flourescent signal from the unknown. The procedure to estimate concentration from flourescence amounts to extending a horizontal line from the observed flourescence and dropping a vertical line at the point of intersection with the standard curve.

For immunoassays, the dose-response curve is monotonically increasing. Interpretations of the 4- and 5PL model parameters are described in Table 1. See Gottschalk and Dunn^{1} for additional details.

Figure 1: Five-parameter logistic dose-response curves and corresponding precision profiles for two different 96-well plates measuring IL7 .

2 The general dose-response model

The general dose-response model is

y_{ij} = f(x_{i}, β) + σg(μ_{i}, θ)ε_{ij}.  (2)

For our purposes, we assume

• y_{ij} is the jth replicate response at the ith dose, j = 1, …, m_{i} and i = 1, …, N where m_{i} is the number of replicates for the ith dose and N is the number of doses.

• x_{i} is the ith dose

• β = (a, b, c, d, g)′

• f(x_{i}, β) = μ_{i} and is the mean response at dose x_{i}. Here f(x_{i}, β) is the function described in Equation 1.

• σ is a an unknown constant

• θ is a vector of parameters for the variance function described by g

• g(μ_{i}, θ) is a variance function that depends on both the mean response at dose i and parameter vector θ

• ε ∼ Normal(0, 1) Model 2 accommodates non-constant variance with

Var(y_{ij}) = σ^{2}g(μ_{i}, θ)^{2}.  (3)

As noted by Gottschalk and Dunn^{1}, the variance at the high-response end of the curve is often three- to four-times larger than at the low end, due to: 1) signal detectors producing noise with a standard deviation proportional to the magnitude of the response; and 2) the kinetics associated with antibody binding are non-linear. For immunoassays, the ‘power of the mean’ variance structure has been shown useful for 4- and 5PL dose-response models^{1}^{,}^{2}, and is modeled as

g(μ_{i}, θ) = μ_{i}^{θ}.  (4)

Here, θ is constrained to be positive.

3 Fitting the standard curve

Typically, the standard curve is fit to ‘background corrected’ response measures. This is achieved by measuring (usually from replicate samples) the flourescent signal from matrix-only samples (i.e. analyte concentration equals 0), and subtracting the average of these readings from all response measures. The standard curve is fit to the backgroundcorrected data using the method of generalized least squares with variance function estimation described in detail by Carroll and Ruppert^{3}, Davidian and Haaland^{4}, and O’Connell, Belanger and Haaland^{2}. Generalized least squares refers to standard weighted least squares with estimated (rather than known) weights. Here, appropriate weights are given by

w_{i} = 1/g^{2}(μ_{i}, θ)  (5)
= 1/μ_{i}^{2θ}.

For completeness, we describe model fitting for the general model with mean function given by Equation 1 and variance function given by Equation 4.

1. Davidian and Haaland^{4} suggest β̂^{(0)} be obtained by minimizing ∑_{i=1}^{N} ∑_{j=1}^{m_{i}}(y_{ij} − μ_{i})^{2}. Note β̂^{(0)} is the least-squares estimate of β. Alternatively, good starting values may already be available for β from the fit obtained from the commercial software. In practice, I have found generating β̂^{(0)} from least-squares estimation to be problematic for 5PL. An ad hoc approach I’ve successfully employed uses the approach outlined below:

(a) a_{0} = max(response)×1.01

(b) d_{0} = min(response)×0.99

(c) g_{0} = 1

(d) Let ỹ = \frac{a−d}{y−d} − 1. Then from Equation 1, log ỹ = b log x − b log c. Let α̂_{0} and α̂_{1} be the fitted intercept and slope from the simple linear regression of log ỹ on log x. Then b_{0} = α̂_{1} and c_{0} = e^{−α̂_{0}/b_{0}}.

2. Obtain θ̂^{(0)} as the slope parameter from a linear regression of log s_{i} on log ȳ_{i}, where s_{i} is the estimated standard deviation of the responses at dose i and ȳ_{i} is the average response at dose i. This is based on the fact that log{Var(y_{ij})}^{1/2} = log σ +θ log μ_{i}.

3. Obtain σ̂^{2(0)} = e^{2α_{0}} where γ_{0} is the estimated intercept from the linear regression in Step 2.

4. Let k = 1.

5. Obtain β̂^{(k)} by minimizing ∑_{i=1}^{N} ∑_{j=1}^{m_{i}} ŵ_{i}^{(k−1)}(y_{ij}−μ_{i})^{2}, where ŵ_{i}^{(k−1)} = 1/g^{2}(μ̂_{i}^{(k−1)}, θ̂^{(k−1)}).

6. From β̂^{(k)}, construct estimates of r_{ij}^{(k)} = y_{ij} − μ̂_{i}^{(k)}.

7. Obtain θ̂^{(k)} by minimizing

∑_{i=1}^{N} ∑_{j=1}^{m_{i}} \{\frac{r_{ij}^{(k)}ˆġ^{(k)}}{g(μ̂_{i}^{(k)}, θ)}\}^{2} ,

where

ˆġ^{(k)} = exp \{\frac{θ}{n} ∑_{i=1}^{N} m_{i} log μ̂_{i}^{(k)}\} .

8. Obtain

σ̂^{2(k)} = \frac{1}{n − p} ∑_{i=1}^{N} ∑_{j=1}^{m_{i}} \frac{(r_{ij}^{(k)})^{2}}{g^{2}(μ̂_{i}^{(k)}, θ̂^{(k)})},

where p is the number of parameters in the mean model (p = 4 or 5 for 4- and 5PL, respectively).

Iterate steps 5 through 8 to convergence. Call the estimates β̂_{GLS}, θ̂_{GLS} and σ̂2_{GLS}. The variance-covariance matrix of β̂_{GLS} is

Σ_{GLS} = σ^{2}(X′G^{−1}X)^{−1}  (6)

where X is an n × p matrix (n = ∑_{i} m_{i}) such that the m_{i} rows corresponding to the ith dose contain the row vector

(\frac{∂}{∂β_{1}} f(x_{i}, β), …, \frac{∂}{∂β_{p}} f(x_{i}, β))′ ,

and G is an n×n diagonal matrix with the diagonal elements of the m_{i} rows corresponding to the ith dose equal to g^{2}(μ_{i}, θ). Equation 6 is estimate by replacing σ^{2} with σ̂2_{GLS}, evaluating X at β̂_{GLS}, and evaluating G at β̂_{GLS} and θ̂_{GLS}.

For the 5PL model, the partial derivatives defining X are as follows:

• ∂f(x_{i}, β)/∂a = (1 + (\frac{x_{i}}{c})^{b})^{−g}

• ∂f(x_{i}, β)/∂b = (a − d)(−g) (1 + (\frac{x_{i}}{c})^{b})^{−g−1} (\frac{x_{i}}{c} )^{b} ln (\frac{x_{i}}{c} )

• ∂f(x_{i}, β)/∂c = (a − d)(−g) (1 + (\frac{x_{i}}{c})^{b})^{−g−1} (\frac{−bx_{i}^{b}}{c^{b+1}} )

• ∂f(x_{i}, β)/∂d = 1 − (1 + (\frac{x_{i}}{c})^{b})^{−g}

• ∂f(x_{i}, β)/∂g = −(a − d) (1 + (\frac{x_{i}}{c})^{b})^{−g} ln (1 + (\frac{x_{i}}{c})^{b})

4 Constructing the precision profile

Background-corrected flourescence readings from subject samples are converted to concentration estimates based on the inverse of Equation 1,

x = f^{−1}(y) = c [(\frac{a − d}{y − d})^{1/g} − 1]^{1/b} .  (7)

Specifically, the estimated concentration for the ith subject’s jth response is x̂_{ij} = f^{−1}(y_{ij}, β̂). If concentration estimates are based on the backfit of the average of replicate flourescence readings (which is typically the case for commercially available software), then a subject-level concentration estimate is x̂_{i} = f^{−1}(ȳ_{i}, β̂), where ȳ_{i} = ∑_{j=1}^{n_{i}} y_{ij} and n_{i} is the number of replicates for the ith subject. Table 4 shows a sample list of (fictitious) concentrations from a multiplex immunoassay. There are three types of ‘data’: the values that look ‘real’ (IDs 2, 3 and 4), the ‘OOR <’ values (IDs 1, 6 and 7), and the starred values (ID 5). ‘OOR’ stands for out of range and indicates flourescence readings smaller than the lower asymptote. Flourescence readings exceeding the upper asymptote are denoted as ‘OOR >’ in the output. In our experience, the majority of out-of-range values occur at lower concentrations. We therefore restrict our discussion to this case, although our results generalize to the analysis of values exceeding the upper asymptote as well. Since it is impossible to obtain a concentration estimate for these response values, the software simply flags them accordingly. The starred values indicate observations that extrapolate beyond the standard range, a practice that is ill-advised.

Correct analysis of subject concentration measures requires determinination of minimum and maximum acceptable concentrations (MinAC and MaxAC, respectively). These thresholds represent concentrations such that values less than (greater than) MinAC (MaxAC) are considered too variable to be reported reliably. Minimum and maximum acceptable concentrations are typically defined as those values for which the coefficient of variation (CV) of the estimated concentration equals a pre-specified cutoff. Figure 1 shows the precision profiles associated with the fitted 5-PL dose-response models for IL7 for two different 96-well plates. The precision profile is a plot of the percent CV in estimated concentration versus concentration, and is typically ‘U-shaped’ for sigmoidal dose-response models. The shape reflects the fact that small changes in flourescent signals in the flat regions of the curve result in highly variable predicted concentrations, while responses in the linear range of the curve predict concentrations with considerably less variability. Traditional thresholds for MinAC and MaxAC are 10% CV (limit of quantitation) and 30% CV (limit of detection) (include references here from the interna-

Table 1: Description of parameters for 4- and 5PL dose-response models.

Parameter 4PL  5PL
a  Asymptote at high dose Asymptote at high dose
d  Asymptote at low dose Asymptote at low dose
b Slope parameter 1: Controls rate of approach to a asymptote

2: bg controls rate of approach to d asymptote

c ED50 Jointly with b and g estimates ED_{ν},

where ED_{ν} = c((100/ν)^{1/g} − 1)^{1/b}

g  1: Controls asymmetry

2: bg controls rate of approach to d asymptote

Table 2: Sample multiplex concentration estimates.

Sample Concentration (pg/ml)
1  OOR <
2  2.65
3  4.6
4  1.09
5  *0.70
6  OOR <
7  OOR <

tional standards literature). However, it has been our observation that thresholding the data at 30% CV or lower can result in little data available for analysis. Gottschalk and Dunn^{5} use a 50% CV threshold in their examples, but make no formal recommendation for establishing MinAC and MaxAC.

Construction of the precision profile requires estimation of Var(x̂). Specifically

Var(x̂) ≐ (\frac{∂f^{−1}}{∂y} )^{2} Var(y) + (\frac{∂f^{−1}}{∂β} )′ Var(β̂) (\frac{∂f^{−1}}{∂β} ) ,  (8)

where • Var(y) is given by Equation 3.

• Var(β̂) is given by Equation 6.

• \frac{∂f^{−1}}{∂y} = \frac{c}{b} [(\frac{a−d}{y−d})^{1/g} − 1]^{\frac{1}{b} −1} × \frac{1}{g} (\frac{a−d}{y−d})^{\frac{1}{g} −1} × \frac{d−a}{(y−d)^{2}}

• \frac{∂f^{−1}}{∂β} = (\frac{∂f^{−1}}{∂a} , \frac{∂f^{−1}}{∂b} , \frac{∂f^{−1}}{∂c} , \frac{∂f^{−1}}{∂d} , \frac{∂f^{−1}}{∂g} )′ with

▷ \frac{∂f^{−1}}{∂a} = \frac{c}{b} [(\frac{a−d}{y−d})^{1/g} − 1]^{\frac{1}{b} −1} × \frac{1}{g} (\frac{a−d}{y−d})^{\frac{1}{g} −1} × (\frac{1}{y−d})

▷ \frac{∂f^{−1}}{∂b} = −c [(\frac{a−d}{y−d})^{1/g} − 1]^{\frac{1}{b}} × log [(\frac{a−d}{y−d})^{\frac{1}{g}} − 1] × b^{−2}

▷ \frac{∂f^{−1}}{∂c} = [(\frac{a−d}{y−d})^{1/g} − 1]^{\frac{1}{b}}

▷ \frac{∂f^{−1}}{∂d} = \frac{c}{b} [(\frac{a−d}{y−d})^{1/g} − 1]^{\frac{1}{b} −1} × \frac{1}{g} (\frac{a−d}{y−d})^{\frac{1}{g} −1} × \frac{a−y}{(y−d)^{2}}

▷ \frac{∂f^{−1}}{∂g} = \frac{c}{b} [(\frac{a−d}{y−d})^{1/g} − 1]^{\frac{1}{b} −1} × [− (\frac{a−d}{y−d})^{\frac{1}{g}} ] × log (\frac{a−d}{y−d}) × g^{−2}.

Estimation of Equation 8 is achieved by substituting σ̂^{2}g^{2}(x, θ̂) for Var(y) and evaluating Equation 8 for a grid of x̂ and y values covering the domain of the standards. The estimated percent coefficient of variation is \sqrt{V̂ar(x̂)}/x̂ × 100 and is plotted against x̂ for the entire domain of standard concentrations as the precision profile. Precision profiles for two different standard curves are shown in Figure 1. A threshold of 50% CV is drawn indicating MinAC and MaxAC values.

5 Simulation study

The purpose of the current study is to investigate operating characteristics of various choices of MinAC and MaxAC with the goal of establishing guidelines for investigators implementing this technology. We conducted the following simulation study.

1. Assume the uncorrected dose-response model is

z_{ij} = f(x_{i}, β) + k + σg(μ_{i}, θ)ε_{ij},  (9)

where z_{ij} is the uncorrected fluorescence measure corresponding to the jth replicate of the ith concentration. Here, k is the background parameter reflecting a vertical shift in fluorescence readings based on blanks. Fix (β, k, σ, θ)′ = (a, b, c, d, g, k, σ, θ). These are treated as the true dose-response parameters. In (9), f(x_{i}, β) = μ_{i} is given by Equation 1, and g(μ_{i}, θ) = μ_{i}^{θ} and is given by Equation 4.

2. Fix nine standards - eight with known concentrations and one blank (i.e. 0 concentration). Call these standards s_{0}, s_{1}, …, s_{8}, with s_{0} indicating the blank.

3. Generate two blank fluorescent readings based on the following assumed model:

z_{0j} = k + σε_{ij},

where j = 1, 2. Compute the average of these readings, z̄_{0} = (z_{01} + z_{02})/2. This average is our estimate of k; that is to say, k̂ = z̄_{0}.

4. For each of the standards s_{1}, …, s_{8}, use Equation 9 with the true parameter values to obtain duplicate uncorrected flourescent readings, resulting in the sixteen data pairs

((s_{1}, z_{11}), (s_{1}, z_{12}), (s_{2}, z_{21}), (s_{2}, z_{22}), …, (s_{8}, z_{81}), (s_{8}, z_{82})) .

5. Use the background estimate obtained from Step 3 to ‘correct’ the fluorescent readings obtained in Step 4. Specifically, let y_{ij} = z_{ij} −k̂. The corrected readings, {y_{ij}}, will be used to fit the standard curve and to conduct inference.

6. Using the data pairs from Step 5, fit the 5PL standard curve given by Equation 2 using the methods described in Section 3. Generate the estimated dose-response parameters (β̂, σ̂, θ̂)′ = (â, b̂, ĉ, d̂, ĝ, σ̂, θ̂).

7. Construct the precision profile for the standard curve constructed in Step 6 using methods described in Section 4. Identify MinAC and MaxAC corresponding to 10, 20, 30, 40, and 50% CV. Use Equation 1 to find corresponding values of MinAFl and MaxAFl (AFL = acceptable fluorescence). That is to say, MinAFl = f(MinAC, β̂) and MaxAFl = f(MaxAC, β̂).

8. Let c_{i} be the true concentration for the ith subject (i = 1, …, n) where the c_{i} arise

from the model

log(c_{i}) = κ + αg_{i} + ν_{i},

so that

c_{i} = exp{κ + αg_{i} + ν_{i}}.

Here, κ is the overall mean, g_{i} is a group indicator for the ith subject, α is the group-effect, and ν_{i} ∼ Normal(0, σ_{ν}^{2}). If α = 0, the marker is not differentially expressed between the two groups. Generate n c_{i} values.

9. Use Equation 9 with the true parameter values to obtained duplicate uncorrected

subject fluorescent readings for each true concentration. The resulting pairs are

((c_{1}, z_{11}), (c_{1}, z_{12}), …, (c_{n}, z_{n1}), (c_{n}, z_{n2})) ,

where n is the number of subjects.

10. From the {z_{ij}} obtained in Step 9, construct corrected subject fluorescent readings {y_{ij}}, where y_{ij} = z_{ij} − k̂.

11. Using the corrected duplicate subject fluorescent readings, construct a censoring indicator τ_{ij} such that τ_{ij} = 1 if y_{ij} > MinAFl and τ_{ij} = 0 if y_{ij} ≤ MinAFl, where i = 1, …, n and j = 1, 2.

12. Using Equation 7 and β̂ obtained in Step 6, generate observed duplicate concentrations, x_{ij}, as follows:

x_{ij} = \{ \begin{matrix}f^{−1}(y_{ij}, β̂) & if τ_{ij} = 1 \\ MinAC & if τ_{ij} = 0 \end{matrix}

The observed data for analysis is represented by the following ordered pairs: ((x_{11}, τ_{11}), (x_{12}, τ_{12}), …, (x_{n1}, τ_{n1}), (x_{n2}, τ_{n2})) .

13. Using the observed data from Step 12, we assume

log(x_{ij}) = η + γg_{i} + ζ_{i} + δ_{ij}

where • η is the overall mean

• γ is the group effect

• g_{i} is a group indicator

• ζ_{i} is a subject-specific random effect such that ζ_{i} ∼ Normal(0, σ_{ζ}^{2})

• δ_{ij} ∼ Normal(0, σ_{δ}^{2})

• ζ_{i} and δ_{ij} are independent Then the likelihood is given by

L(Ψ|X) = ∏_{i=1}^{n} [∫_{−∞}^{∞} {∏_{j=1}^{2} h(log(x_{ij})|ζ_{i})^{τ_{ij}}H(log(x_{ij})|ζ_{i})^{1−τ_{ij}} r(ζ_{i})} dζ_{i}] ,(10)

where • Ψ = (η, γ, σ_{ζ}^{2}, σ_{δ}^{2})′

• X = (x_{11}, x_{12}, …, x_{n1}, x_{n2})′

• h(log(x_{ij})|ζ_{i}) is a normal probability density function (conditional on the random effect) with mean η + γg_{i} + ζ_{i} and variance σ_{δ}^{2} • H(log(x_{ij})|ζ_{i}) is the normal cumulative density function (conditional on the

random effect) corresponding to h(log(x_{ij})|ζ_{i})

• r(ζ) is the distribution of the random effect and is a normal probability density function with mean 0 and variance σ_{ζ}^{2}

The MLE of Equation 10 is used to conduct group-level inference.

14. Repeating Steps 4 through 13 for the same values of (a, b, c, d, g, k, σ, θ)′ (the true standard curve parameters) and the same values of (κ, α, σ_{ν}^{2})′ (the true subject log concentration parameters) is equivalent to a complete analysis of a single plate. If multiple plates are used in an experiment, then Step 13 should be conducted from the aggregation of subject data from multiple plates.

References

1. P.G. Gottschalk and J.R. Dunn. The five-parameter logistic: A characterization and comparison with the four-parameter logistics. Analytical Biochemistry, 343:54–65, 2005.

2. M.A. O’Connell, B.A. Belanger, and P.D. Haaland. Calibration and assay development using the four-parameter logistic model. Chemometrics and Intelligent Laboratory Systems, 20:97–114, 1993.

3. R.J. Carroll and D. Rupert. Transformation and Weighting in Regression. Chapman & Hall, New York, 1988.

4. M. Davidian and P.D. Haaland. Regression and calibration with nonconstant error variance. Chemometrics and Intelligent Laboratory Systems, 9:231–248, 1990.

5. P.G. Gottschalk and J.R. Dunn. Determining the error of dose estimates and minimum and maximum acceptable concentrations from assays with nonlinear dose-response curves. Computer Methods and Programs in Biomedicine, 80:204–215, 2005.