Prediction of lethal/effective concentration/dose in the presence of multiple auxiliary covariates and components of variance

Abstract

Predictors of the percentile lethal/effective concentration/dose are commonly used measures of efficacy and toxicity. Typically such quantal‐response predictors (e.g., the exposure required to kill 50% of some population) are estimated from simple bioassays wherein organisms are exposed to a gradient of several concentrations of a single agent. The toxicity of an agent may be influenced by auxiliary covariates, however, and more complicated experimental designs may introduce multiple variance components. Prediction methods lag examples of those cases. A conventional two‐stage approach consists of multiple bivariate predictions of, say, medial lethal concentration followed by regression of those predictions on the auxiliary covariates. We propose a more effective and parsimonious class of generalized nonlinear mixed‐effects models for prediction of lethal/effective dose/concentration from auxiliary covariates. We demonstrate examples using data from a study regarding the effects of pH and additions of variable quantities 2',5'‐dichloro‐4'‐nitrosalicylanilide (niclosamide) on the toxicity of 3‐trifluoromethyl‐4‐nitrophenol to larval sea lamprey (Petromyzon marinus). The new models yielded unbiased predictions and root‐mean‐squared errors (RMSEs) of prediction for the exposure required to kill 50 and 99.9% of some population that were 29 to 82% smaller, respectively, than those from the conventional two‐stage procedure. The model class is flexible and easily implemented using commonly available software.

Quantal‐response bioassay, Generalized nonlinear mixed model, Marginal prediction, 3‐Trifluoromethyl‐4‐nitrophenol, Sea lamprey

INTRODUCTION

Quantities such as the percentile lethal/effective dose/concentration predict the exposure required to evoke some quantal response (e.g., death) in 100p*% of some population, where 0 < p* < 1. Such prevalence‐predicting exposures, herein denoted LC100p* for simplicity, are commonly used end points in toxicity studies and regulatory assessments. The median lethal concentration (LC50), for example, is defined by p* = 0.5. Values of LC100p* are not observable as primary data. Rather, they must be predicted from the inversion of models that predict the probability of the quantal response from exposure dose or concentration. Typically, LC100p* is predicted from relatively simple bivariate bioassays wherein organisms are exposed, often groupwise, to one of several concentrations of a single agent, and much attention has been devoted to that prediction [1–6].

The toxicity or efficacy of an agent also might be affected by multiple factors, including the presence of bivalent cations [7], pH [8], or dissolved organic carbon [9]. In such cases, it is desirable to predict LC100p* over the ranges of those auxiliary covariates, but the development of methods has lagged behind need. Conventionally, an intuitively reasonable two‐stage approach has been used wherein multiple LC100p* values—each from a specific level of an ancillary covariate—are predicted using conventional bivariate models and subsequently treated as observed data in second‐stage regression models that include the auxiliary covariates [9–11]. That conventional two‐stage approach lacks parsimony and coherence, however, because estimation of LC100p* and prediction over the auxiliary covariates are separated. Furthermore, the uncertainty of the LC100p* estimates may be misrepresented if the designs of the studies entail variance components that cannot be accommodated using that conventional two‐stage approach. We propose a simple and coherent class of models to overcome those limitations.

The present study was motivated by the need to provide improved estimates of the effects of 2',5'‐dichloro‐4'‐nitrosalicylanilide (niclosamide) on the toxicity of 3‐trifluoromethyl‐4‐nitrophenol (TFM) to larval sea lamprey (Petromyzon marinus) over a wide range of ambient pH. Adult sea lamprey are parasitic on other fishes, and repeated attacks can be lethal. The Laurentian Great Lakes ecosystem was radically altered by high mortality of large‐bodied native fishes following the invasion of sea lamprey after the opening of the Welland Canal, and TFM has been used to control sea lamprey in the Great Lakes of North America since approximately 1959 [12].

The uptake of TFM decreases with increasing pH, because the ionized form of the TFM molecule, which becomes increasingly dominant as ambient pH increases beyond its dissociation constant (pK_a) of 6.07, transfers less readily across the gill membrane than the unionized form [13]. This sensitivity to pH is critically important, because the pH of spawning tributaries varies spatially and temporally. Additions of small amounts of niclosamide are used to boost toxicity, especially in high‐pH streams. Therefore, successful treatment with these agents requires prediction of toxicant formulation and concentration, based on pH, for each application to a spawning tributary.

The treatment of sea lamprey spawning areas is expensive and can have untoward effects. Although TFM and niclosamide are relatively selective for larval sea lamprey, juvenile lake sturgeon (Acipenser fulvescens) also are vulnerable [14]. Therefore, current stream‐treatment protocols require application of the minimum formulation needed to kill 99.9% of the larval sea lamprey (LC99.9) while sparing juvenile lake sturgeon. Accuracy and precision are so important that streamside bioassays have been used to estimate toxicant requirements. The standard practice, however, is to predict treatment requirements based on the results of laboratory toxicity tests.

Conventionally, laboratory estimates of LC100p* have been treated as observed data and regressed on niclosamide concentration and either pH or alkalinity to produce tables of pH/alkalinity and niclosamide‐specific treatment concentrations [10]. Treatment managers measure pH or alkalinity in situ and then use the standard tables to identify suitable toxicant formulations and concentrations.

We propose a parsimonious class of models for prediction of LC100p* and associated measures of uncertainty from complex studies that include multiple auxiliary variables and variance components. The model class is broad and, therefore, requires presentation using a rather high degree of mathematical abstraction. However, we provide SAS^® (SAS Institute, Cary, NC, USA) code (see Appendix [http://dx.doi.org/10.1897/06–630.S1]) to demonstrate some examples.

We demonstrate models using data from a study regarding the toxicity of TFM and niclosamide combinations on larval sea lamprey at multiple values of pH. Estimates of LC100p* may vary among models, and we demonstrate inclusion of that model uncertainty using model averaging. Measures of uncertainty are particularly important for that application, because historically, stream‐treatment managers have made small, ad hoc adjustments to formulations based on intuition and experience. Prediction intervals provide data‐based limits on the range of those ad hoc adjustments.

MATERIALS AND METHODS

Review of bivariate quantal‐response models for estimation of LC100p*

A typical quantal‐response bioassay test consists of some number I of experimental units (EUs), each consisting of n_i ⩾ 1 organisms, that are exposed to an increasing gradient of some agent. Let X_i represent the concentration/dose of the agent to which EU i was exposed, such that X_i ⩾ 0 for all i, and let Y_i represent the number of test organisms out of each n_i that exhibited some quantal response (e.g., death) during the test. Control EUs have X_i = 0. Under the condition that the Y_i represent the sums of independent homogeneous Bernoulli random variables, each follows a binomial distribution having parameter p_i, which represents the proportion of organisms that exhibit the quantal response in the population that produced EU i. Therefore, LC100p* is traditionally predicted from binomial models of the form

where T⁻¹(·) is the inverse of some tolerance function T(·), and β₀ and β₁ are parameters describing the relation between p_i and h_Λ(X_i). We demonstrate the commonly used cumulative standard Gaussian, logistic, and extreme‐value tolerance functions (Table 1). The function h_Λ(·) is a device that expands the range of shapes that T(·) can assume, and is equivalent to transformation of exposure X. Prior data transformation is unnecessary, and instead, we incorporate the Box‐Cox power function

to mimic the effect of a broad class of transformations on T(·). Note that h_Λ(X) = X for Λ = 1. The single component of variance is Var(Y_i) = n_ip_i(1 — p_i). Conventional predictions of LC100p*, denoted X̂_CONV(p*), are obtained by estimating β₀ and β₁, then setting p_i = p* and solving Equation 1 for X to predict the concentration/dose that is expected to produce the quantal‐response rate of p* and given by

For the Gaussian tolerance function and Λ = 1, for example, the conventional estimates of LC50 are then given by

These are essentially the Litchfield‐Wilcoxon type estimates produced by specialized software such as TOXSTAT^® [15] and ToxCalc^<35> [16].

Table 1.

Commonly used tolerance functions T(w) = p and their inverses, T⁻¹(p) = w, used in quantal response bioassays^a

Tolerance function/inverse	T(w)	T−1(p)

^a Here, p and w represent the probability of some quantal response and exposure, respectively. These functions can be generalized by defining w = β₀ + β₁h_Λ(x), where h_Λ(x) is any monotonic transformation function in x.

Conventional two‐stage prediction of LC100p* from auxiliary covariates

We consider the more general case where the efficacy or toxicity of an agent may be affected by one or more auxiliary covariates that are not intrinsically toxic or are confounded with the primary agent by design. Typically, the data are generated from multiple clusters, which are aggregations of EUs that might reasonably be correlated. For example, tests—separated in space or time—might be treated as clusters. The auxiliary covariates may vary among clusters. Let j index the J clusters. Conventional ad hoc prediction of LC100p* from the auxiliary covariates is done in two stages. The first stage requires the fitting of J independent conventional quantal‐response prevalence models, given by

to produce J estimates of LC100p*, denoted by X̂_CONVj(p*). Note that those first‐stage models require estimation of a total of 2J parameters. For the second stage, the J estimates of LC100p* are treated as observed data in generalized linear models of the form

where g²(·) is some variance function [17] defining any heteroscedasticity in Var[X̂_CONVj(p*)] = σ²g²(·) and A_ij(θ₀ + z'_ijθ) is a function of parameter θ₀, a vector of M ancillary covariates z_ij = [Z_ij, …, Z_Mj]' and parameter‐vector θ, where z' denotes the transpose of z. For example, if g²(·) = 1 and A(·) is the identity function, then Equation 5b reduces to the ordinary linear regression model:

where Var[X̂_j(p*)] = Var(ε_j) = σ². The conventional two‐stage predictor of LC100p*, denoted X̂_2STG(p*), is obtained from the fitting of Equation 5b as

Note that the X̂_CONVj(p*) are predicted with error from Equation 5a; therefore, Equation 5b is an errors‐in‐response model. The parameters of such models are not universally identifiable without additional information, but the necessary conditions are well known [18] and satisfied by ordinary linear regression. The second‐stage model requires the estimation of the M + 1 structural parameters in A(·), σ², plus any parameters in g²(·). Therefore, the conventional two‐stage prediction of LC100p* from ancillary covariates requires the estimation of a total of at least 2J + M + 2 parameters.

Table 2.

Formulas for the marginal expectation of LC100p*, denoted X̂_MARG (p*), obtained from solutions to Equation 9 for h₁(X_ij) = X_ij and h₀(X_ij) = log(X_ij + 1)^a

Tolerance function	h1(Xij) = Xij	h0(Xij) = log(Xij + 1)

^a For a single random effect, u_j = u_j, v'_ij∑v_ij = σ².

Careful scrutiny also reveals that this superficially intuitive two‐stage model does not seem to represent any coherent and natural underlying data‐generating process. The first stage mimics various binomial data‐generating processes that are explicitly independent of any auxiliary covariates. In contrast, the second stage represents data generation by auxiliary covariates, thereby contradicting the first‐stage processes. The two‐stage model might yield satisfactory results in some cases, but it apparently lacks a logical theoretical basis.

New models for prediction of LC100p* from ancillary covariates

We propose a general and more parsimonious class of normal‐binomial models given by

where i indexes EUs, j indexes clusters, u_j is a q × 1 vector of random effects for the jth cluster with design vector v_ij; and ∑ denotes the variance‐covariance matrix of the u_j. The u_j are estimated from the data. Here, θ₀ is excluded from the argument of A_ij(·), where it would be confounded with β₀. These models contain multiple components of variance; the conditional variance within clusters is Var(Y_ij|u_j) = n_ijp_ij (1 — p_ij) and the among‐cluster variances and covariances are contained in ∑. Furthermore, they require estimation of only M + 2 + rank(∑) parameters, which generally will be far fewer than the conventional two‐stage approach.

The model given by Equation 7 is conditional on the u_j and produces cluster‐specific predictions P̂_ij. The resulting conditional predictor of LC100p* is given by

Note that if û_j = 0 for all j, then the argument of h⁻¹_Λ(·) in Equation 8 simplifies to the customary estimate of LC100p* (Eqn. 3) multiplied by A_ij(z'_ijθ). The conditional model is the fundamental representation of the data [19], but the resulting cluster‐specific predictions are not useful for prediction of stream‐treatment requirements. For that purpose, the marginal (integrated over all of the cluster‐specific random effects) expectation is needed for prediction of LC100p* outside the specific set of tests and is given by

where F(û) is the multivariate normal cumulative distribution function of û. Algebraic solutions to Equation 9 exist for the most common choices of h₁(X) = X and h₀(X) = log(X + 1) (Table 2). Otherwise, the marginal expected LC100p* must be obtained using numerical integration.

The model defined by Equation 7 is general and can be adapted to accommodate complex study designs. For example, variance might be partitioned across three or more hierarchical levels (e.g., individuals within tanks, among tanks within tests, and among tests) [20] that requires only modification of the specification of the random‐effects formulation v'_iju_j and covariance matrix ∑. Other tolerance functions (e.g., cumulative Weibull) also might be incorporated but are not considered further here.

Alternative formulations

Our model includes random effects on the scale of β₀. Random effects may enter at the scale of other parameters (i.e., β₁). We did not consider those alternatives further for lack of need in our data. The parameters β₀ and β₁ tend to be highly correlated in tolerance functions; therefore, there may be little need or ability to model random effects for both parameters.

Note, from Equation 7, the multiplicative effect of the auxiliary variables z_ij through the product X_ij × [A_ij(z'_ijθ)]⁻¹. Therefore, the auxiliary variables do not exert any identifiable effect on p_ij when X_ij = 0. That formulation is sensible for auxiliary variables that affect the potency of the primary agent X but that otherwise are inert within the range of the data or are confounded with X by design. In contrast, an additive conditional inverse tolerance function, defined by T⁻¹ (p_ij) = β₀ + v'_iju_j + β₁X_ij + A(z'_ijθ), permits the auxiliary variables to exert effects independently of X_ij so that p_ij may vary systematically even when X_ij  0. An additive or hybrid formulation may be an appropriate alternative for studies of multiple toxins. For example, suppose that the concentrations of two toxins X_1ij and X_2ij, each affected by nontoxic covariates z_ij, are crossed in a factorial design. A plausible hybrid formulation for quantal‐response function is

which yields two nonunique marginal LC100p* estimates given by

This formulation is easily extended to the case where X₁ and X₂ are affected by different covariates. We did not consider such additive or hybrid formulations for our application, because TFM and niclosamide concentrations were confounded by the design and pH is nontoxic within the range of the data.

Identification of a form for the effect of the ancillary covariates

An appropriate form for A_ij(z'_ijθ) may be identified from existing knowledge or the relation between the conventional LC50 estimates and the Z_mij in the data. Note that from Equation 8, h_Λ[X̂_j(p* | z_ij, u_j)] is proportional to A_ij(z'_ijθ). The choice p* = 0.5 is convenient for preliminary identification of A_ij(z'_ijθ), but other choices also are viable. The X̂_CONVj(0.5), computed independently from each test using Equation 4, can be viewed as näive preliminary estimates of the X̂_j(0.5 | z_ij,u_j) in Equation 8 to identify a plausible form for A_ij(z'_ijθ). For example where h₁(X) = X and plots of the log[X̂_CONVj(0.5)] versus each of the M ancillary covariates Z_mij reveal linear scatters, then

Careful exploratory analyses will be needed to identify an appropriate form for A_ij(z'_ijθ) in each study.

Parameter estimation

The model given by Equation 7 is an example of a generalized nonlinear mixed‐effects model and can be fitted using maximum pseudolikelihood estimation [21], Bayesian Markov‐Chain Monte Carlo sampling [22], or maximum likelihood estimation [23] using commonly available software. Maximum likelihood estimates can be obtained using SAS NLMIXED software or J.K. Lindsey's repeated package for R [24] (opensource software: http://cran.r‐project.org/ and http://popgen.unimaas.nl/∼jlindsey/rcode.html). We demonstrate model fitting using the SAS NLMIXED package, which deploys the method of Pinheiro and Bates using Gaussian quadrature to integrate the likelihood over the random effects. Standard errors and Wald‐type frequentist prediction intervals on the marginal LC100p* estimates Equation 9 are approximated by NLMIXED using the delta method [25] and prediction standard errors [26].

Acute toxicity tests on sea lamprey

The motivating example for the present study was a set of 38 tests designed to assess the effects of pH and niclosamide augmentation on the toxicity of TFM. The test agents were isopropanol solution of the sodium salt formulation of TFM, with an average concentration of 38% active ingredient (Lamrecid^®; Farbwerke Hoechst AG, Frankfurt, Germany), and the aminoethanol salt formulation of niclosamide, with an average niclosamide concentration of 59% (Bayluscide^®; ProServe, Menphis, TN, USA). Sea lamprey ammocetes (80–120 mm) were collected from the Ford River (MI, USA) and held for a minimum of 14 d before testing to eliminate any effect of delayed handling mortality from the toxicity trials.

All toxicity testing followed standard protocols [27,28]. Thirty‐eight tests were performed as 19 paired runs using a single, continuous‐flow, serially diluting delivery system similar to that described by Garton [29]. The delivery system consisted of two parallel but independently controllable series, each consisting of 10 glass exposure tanks. Tanks were the EUs and, in each series, were independently randomized to treatments. Each test consisted of gradients of TFM concentration sufficient to produce mortalities from 0 to 100% at each pH (7, 8, and 9) (Table 3). For each run, one series was randomly selected for augmentation with niclosamide, and the other served as a niclosamide‐free test. The niclosamide augmentations ranged from 0.125 to 2.000% of the TFM mass concentration. Therefore, TFM concentration and niclosamide addition were multiplicatively confounded in this design. Ten healthy sea lamprey ammocetes were randomized to each tank 24 h before the start of each test. The 19 matched pairs of tests were performed in a random temporal sequence. Each test consisted of 12 h of exposure, followed by 12 h of postexposure observation for delayed mortality. Dead sea lamprey were removed hourly from each tank.

Well water was used for all exposures. Water temperature (mean, 12.5°C; range, 11.9–13.4°C), pH, and dissolved O₂ concentration (mean, 9.1 mg/L; range, 8.4–9.6 mg/L) were measured hourly in the tanks during the exposures. Total alkalinity (mean, 108 mg/L as CaCO₃; range, 82–129 mg/L of CaCO₃) and total hardness (mean, 150 mg/L as CaCO₃; range, 145–159 mg/L as CaCO₃) were measured every 4 h in the control tanks [30]. The desired pH of the test water was obtained by metering dilute acid (0.5 M HCl) or base (0.5 M NaOH) into the diluter headbox using a Masterflex^® peristaltic pump (Cole‐Parmer, Vernon Hills, IL, USA). Water samples were collected hourly, and TFM concentrations were measured using a Beckman DU‐640 spectrophotometer (Beckman Coulter, Fullerton, CA, USA). Niclosamide concentrations were measured every 2 h by high‐performance liquid chromatography using a Millennium Series 2010 LC Module I Plus (Waters, Milford, MA, USA).

Table 3.

Realized 3‐trifluoromethyl‐4‐nitrophenol (TFM) concentrations to which sea lamprey (Petromyzon marinus) ammocetes were exposed^a

	TFM (mg/L)
Experimental unit	Gradient 1 (pH 7)	Gradient 2 (pH 8)	Gradient 3 (pH 9)
1	0.00	0.00	0.00
2	0.09	0.19	0.26
3	0.37	0.26	0.48
4	0.41	0.54	0.75
5	0.80	0.82	1.11
6	1.06	1.02	2.25
7	1.35	1.70
8	1.88	2.15

^a Tests were defined as combinations of pH and niclosamide treatment. Each test required a series of six or eight tanks comprising pH‐specific TFM gradients. At each pH, paired tests were conducted in which niclosamide was added to each tank of one series, and the other series served as a niclosamide‐free treatment. The target niclosamide additions were 0.125, 0.250, 0.500, 0.750, 1.000, 1.500, and 2.000% of the mass of the TFM concentration; however, the 0.125 and 1.500% niclosamide additions were omitted at pH 9 for lack of ammocetes. Therefore, 19 combinations of pH and niclosamide augmentation plus 19 niclosamide‐free tests were used. Actual niclosamide additions varied somewhat from the nominal targets and were used in all analyses. The numbering of experimental units (tanks) is arbitrary; treatments were randomized to tanks in each test.

Application of the new models to the sea lamprey toxicity tests

Each test was sufficient to independently compute X̂_CONVj(p*) for each particular combination of pH and niclosamide augmentation. Preliminary analyses of the 38 X̂_CONVj(0.5) indicated that the LC50 of TFM varied approximately exponentially with a quadratic polynomial in niclosamide augmentation Z₁ and pH Z₂ so that

was a reasonable form. Combining Equations 7 and 10, appropriate models that include a single random effect (q = 1) are given by

where σ² = ∑ denotes among‐cluster variation in β₀. From Equation 9, the formulas for X̂_MARG(p*)for our choices of T(·) and h_Λ(X_ij) are given by substitution of the right‐hand side of Equation 10 for A(z'_ijθ) in Table 2. For computational convenience, we centered pH at 8 so that Z₂ = pH — 8. In preliminary analyses, we fitted the nonlinear mixed‐effects model using tests, runs (matched pairs of tests), and run × tank combinations as clusters. Based on Akaike's information criterion (AIC; smaller values are better) [31], the choice of tests for cluster definition produced, by far, the best models, and we adopted that choice for subsequent analyses. We fitted logistic, probit, and cumulative extreme‐value model variations for h₁(X_ij) = X_ij and h₀(X_ij) = log(X_ij + 1). Example SAS NLMIXED code for fitting the models given by Equation 11 and marginal prediction of LC100p* is given in the online‐only Appendix (http://dx.doi.org/10.1897/06–630.S1).

In practice, uncertainty about model choice cannot be eliminated. We accommodated model uncertainty by frequentist AIC‐weighted averaging [32,33] of the resulting six marginal predictors X̂_MARG(p*)and their 95% prediction limits.

Model comparison

For comparison we mimicked the conventional two‐stage practice used in the sea lamprey control program. We first obtained the conventional estimates of LC50 and LC99.9 from Equation 5a. For the second stage, we fitted the generalized linear model given by

to obtain X̂_2STG(p*). Each of the 38 X̂_CONVj(0.5) required estimation of two parameters, and the secondary model based on Equation 10 required estimation of six additional parameters (θ₀, …, θ₄, and σ²). The resulting 39 models can be viewed as a single, 82‐parameter composite model. We obtained the AIC of the 82 parameter model as the sum of the AIC scores of the 38 fittings of Equation 1 plus the AIC from the nonlinear regression of LC50 estimates on pH and niclosamide augmentation. We compared the performance of that 82 parameter composite model with our seven‐parameter, nonlinear mixed‐effects probit model using AIC; differences larger than 10 may be interpreted as strong evidence [33] of the superiority of our nonlinear mixed‐effects model.

We compared the predictive performance of X̂_MARG(p*)and X̂_2STG(p*) using Monte Carlo simulation. We defined a fully determined hypothetical population based on the predictor variables and parameter estimates from the sea lamprey study. Therefore, the true population LC100p* values for each cluster, denoted X_j(p*), were fixed and known in the simulation. We used Equation 11 with Λ = 1 and the Gaussian tolerance function as our mortality‐generating model. For each of σ² = 0 and σ² = 1.24, we sampled 1,000 replicate sets of mortality counts Y_ij, i = 1, …, 10, and j = 1, …, 38, using independent normal and binomial random‐number generators for u_ij and Y_ij, respectively. The simulation in which σ² was fixed at zero enables comparison without risk of misspecification of amongcluster random variation. We predicted LC50, LC99.9, and their Wald‐type prediction intervals for each cluster within each replicate using X̂_MARG(p*) and X̂_2STG(p*). We estimated the predictive accuracy of X̂_MARG(p*) and X̂_2STG(p*) using Bias_E = (1/J) ∑^J_j=1 [X̂_Ej(p*) — X_j(p*)], where E ε {MARG, 2STG}, and we estimated their precision using RMSE of prediction (RMSE_E = for each Monte Carlo replicate. We estimated coverage probabilities of the 95% Wald‐type prediction intervals as Pr(coverage) = (1/1,000J) ∑^1,000_r=1 ∑^J_j=1  I[X_jr(p*) ε {S_Ejr}] where I(·) ε {0, 1} is an indicator function for its argument and {S_Ejr} is the estimated Wald‐type prediction interval for X̂_Ejr(p*). Perfect prediction intervals have coverage probabilities that achieve their nominal values.

Table 4.

Parameter estimates and Akaike's information criterion (AIC) scores for the six nonlinear mixed‐effects models of effects of 3‐trifluoromethyl‐4‐nitrophenol (TFM) concentration, niclosamide fraction, and pH on mortality of juvenile sea lampreys (Petromyzon marinus)^a

	hΛ(X) = [TFM]			hΛ(X) = log([TFM] + 1)
Parameter	Logistic	Gaussian	CEV	Logistic	Gaussian	CEV
β0	‐18.53 (1.12)	‐10.13 (0.56)	‐12.75 (0.71)	‐29.29 (1.76)	‐15.66 (0.83)	‐19.27 (1.08)
β1	9.55 (0.62)	5.22 (0.31)	6.29 (0.38)	27.25 (1.68)	14.61 (0.81)	17.47 (1.03)
θ1	−0.85 (0.09)	−0.84 (0.09)	−0.83 (0.09)	−0.56 (0.06)	−0.55 (0.06)	−0.55 (0.06)
θ2	1.30 (0.02)	1.30 (0.02)	1.30 (0.02)	0.79 (0.01)	0.79 (0.02)	0.78 (0.02)
θ3	0.22 (0.04)	0.22 (0.04)	0.22 (0.04)	0.16 (0.03)	0.15 (0.03)	0.15 (0.03)
θ4	0.18 (0.04)	0.18 (0.04)	0.18 (0.04)	−0.03 (0.02)	−0.03 (0.02)	−0.02 (0.02)
σ2	4.12 (1.16)	1.24 (0.34)	2.01 (0.55)	4.20 (1.21)	1.25 (0.36)	1.97 (0.56)
AIC	332.0	332.6	362.4	347.9	354.4	374.2

^a The models are the combinations of two functions of TFM concentration and the logistic, Gaussian, and cumulative extreme‐value (CEV) tolerance functions. Values in parentheses represent the standard error.

RESULTS AND DISCUSSION

Partial kills of larval sea lamprey were observed in fewer than 20% of the EUs. The remaining EUs were nearly equally split between 0 and 100% mortality rates. The large proportions of data at the survival extremes aided precision and the ability to estimate LC100p* for extreme values of p*. These data, however, also challenge estimation of the variance of X̂_CONVj(p*). For example, Fieller's theorem, which is the commonly used basis for estimation of confidence intervals for X̂_CONVj(p*), lacked real‐valued solutions in 16 of the 38 tests and, therefore, was not useful for these data. The delta‐method approximation also failed to produce confidence intervals for X̂_CONVj(p*) from 8 of the 38 tests, because positive‐definite estimates of covariance matrices could not be obtained.

Predictions of LC100p* of TFM for sea lamprey

Our nonlinear mixed‐effects models based on TFM concentration fit better than models based on log(TFM), and the logistic and Gaussian tolerance functions fit better than the cumulative extreme‐value function (Table 4). The AIC‐weighted marginal predictor X̂_MARG(p*) was dominated by the logistic and Gaussian tolerance functions, because those were vastly superior for our data based on AIC. All seven parameters were important and estimated with reasonable precision. The among‐cluster variances, σ², ranged from 1.2 to 4.2, and their confidence intervals never included zero, confirming the importance of modeling that source of variation. No evidence of any systematic departure from the data was found

Although pH in the range of 7 to 9 is not toxic to larval sea lamprey, pH exerts a strong influence on the toxicity of TFM. Niclosamide augmentation, within the range used in the present study, has relatively little effect on the toxicity of TFM to larval sea lamprey at pH 7, but such augmentation substantially enhances toxicity at pH values greater than 8 (Fig. 1). Generally, little difference was found between LC50 and LC99.9 and, hence, little margin for error in the application of those formulations. These results underscore the importance of precise in situ measurement of pH before each treatment of sea lamprey spawning tributaries to identify the most efficient and effective combination of TFM and niclosamide. Estimates of precision are critical to both regulatory requirements [8] and identification of the plausible range of outcomes from any treatment application. The widths of the pointwise 95% intervals on our marginal predictions of LC50 and LC99.9 were small relative to those predictions (Fig. 1). Our model provides precise predictions of LC50 and LC99.9 at any pH from 7 through 9 and any addition of niclosamide from 0 through 2% of the mass of TFM.

Comparison of the new and conventional predictors

The two‐stage predictor based on the Gaussian tolerance function X̂_2STG(0.5) required estimation of 82 parameters and had an AIC score of 373.0. Based on the resulting AIC difference of 40.4, X̂_MARG(p*) for the Gaussian tolerance function (AIC = 332.6) (Table 4) outperformed X̂_2STG(p*) by a wide margin. The superior AIC of the generalized nonlinear mixed model was entirely the result of its parsimony. The maximized log‐likelihood of the two‐stage model was −104.5 and exceeded the value of −186.5 from the generalized nonlinear mixed model. It does not matter, however, whether an AIC difference between two models is caused by the number of parameters or by the values of the maximized log‐likelihoods, because AIC is an omnibus measure that approximates the Kullback‐Liebler distance between an approximating model and the unknown truth it represents. We should expect better predictions from those models that have smaller AIC values.

The marginal predictor based on the generalized nonlinear mixed probit model, X̂_MARG(p*), clearly outperformed the conventional two‐stage predictor X̂_2STG(0.5) in the Monte Carlo simulations. The 95% prediction intervals for X̂_2STG(p*) were too narrow (Table 5). In contrast, the 95% intervals for X̂_MARG(p*) achieved the nominal coverage for σ² = 0 but were somewhat too narrow for σ² = 1.24. The median RMSEs of prediction for X̂_MARG(p*) were 43 and 83% smaller than those for X̂_2STG(p*) for p* = 0.5 and p* = 0.999, respectively, when σ² = 0 (Fig. 2). The marginal predictor X̂_MARG(p*) was unbiased, whereas X̂_2STG(0.999) underestimated the LC99.9 of the TFM‐niclosamide formulations by approximately 0.26 mg/L. The RMSEs of prediction forX̂_2STG(p*) were more variable and, sometimes, unacceptably large at σ² = 1.24 (Fig. 3). Again, the median RMSEs of prediction were 29 and 45% smaller for p* = 0.5 and p* = 0.999, respectively, for X̂_MARG(p*) than for X̂_2STG(p*). The biases of both X̂_2STG(p*) and X̂_MARG(p*) tended to be negligibly small, but X̂_2STG(0.999) usually underestimated LC99.9. Furthermore, the small RMSE and negligible bias of X̂_MARG(0.999) in the Monte Carlo simulation for σ² = 1.24 clearly indicate that LC99.9 can be predicted with good accuracy and precision from trials similar to the sea lamprey study.

Fig. 1.

Predictions of concentrations lethal to 50 and 99.9% of exposed organisms (LC50 and LC99.9, respectively) of 3‐trifluoromethyl‐4‐nitrophenol (TFM) for larval sea lamprey (Petromyzon marinus) based on Akaike's information criterion (AIC)‐weighted averages of marginal predictors X̂_MARG(p*) and the widths of their pointwise 95% prediction intervals (PIs; 95% PI width). The model averaging spanned the six combinations of the logistic, Gaussian, and cumulative extreme‐value tolerance functions each based on [TFM] and log([TFM] + 1).

Table 5.

Coverage probabilities for 95% prediction intervals on the concentrations expected to kill 50 and 99.9% of the organisms (LC50 and LC99.9, respectively) in arbitrary populations based on X̂_2STG(p*) and X̂_MARG(p*)^a

		Predictor
Among‐cluster variance	Predicted exposure	X̂2STG(p*)	X̂MARG(p*)
σ2 = 0	LC50 (p* = 0.5)	0.88	0.95
	LC99.9 (p* = 0.999)	0.74	0.94
σ2 = 1.24	LC50 (p* = 0.5)	0.88	0.90
	LC99.9 (p* = 0.999)	0.84	0.90

^a Coverage probabilities were estimated from Monte Carlo simulations consisting of 1,000 random replicates of the sea lamprey (Petromyzon marinus) trials for each choice of σ². The correct coverage probability for a 95% prediction interval is 0.95. Nominal 95% prediction intervals having coverage probability less than 0.95 err by being too narrow.

Additional considerations

We relied on the asymptotic properties of the delta‐method approximation of 95% prediction intervals. Based on our Monte Carlo simulation, that approximate method performed extremely well for X̂_MARG(p*) when σ² = 0, but it was unlikely to achieve the nominal coverage for the sea lamprey data for which σ = 1.24. The delta‐method approximation also performed reasonably well for estimation of confidence intervals for the median lethal dose in simulation studies for which J ⩾ 30 [34] and may be adequate for studies that include relatively large numbers of EUs. Clearly, however, room exists for improvement, and smaller studies may require other approaches [35,36]. A further refinement or evaluation of prediction intervals was beyond the scope of the present study.

Bias also is a critical issue. We examined estimation bias, assuming correct model specification, in our Monte Carlo simulations, and X̂_MARG(p*) performed well. Bias also can result from model misspecification, particularly the form for A(z'θ) in Equations 5b and 7. We explored many parametric alternatives for A(z'θ) that did not perform as well as our final choice. Those alternatives included variations of the simple linear form A_ij (z'θ) = θ₁Z_1ij + θ₂Z_2ij + θ₃Z_1ijZ_2ij, exponential functions of those forms, and exponentials of a broad set of fractional polynomials [37]. The P̂_ij produced by our choice for A(z'θ) showed no evidence of systematic bias in our data. The choice of A(z'θ), however, merits careful investigation in each study, and the value of careful exploratory analyses cannot be understated. Room exists here for additional development. For example, incorporation of nonparametric smoothers [38] for estimation of A(z) might more automatically minimize the risk of misspecification bias.

Fig. 2.

Performance based on precision (root‐mean‐squared error [RMSE]) and accuracy (Bias) of the conventional two‐stage (X̂_2STG) and the marginal mean (X̂_MARG) predictors of concentrations lethal to 50 and 99.9% of exposed organisms (LC50 and LC99.9, respectively). Performance was estimated from 1,000 Monte Carlo replications of the 38 clusters of the sea lamprey (Petromyzon marinus) trials for which the among‐cluster variance was fixed at σ² = 0. Boxes enclose the 25th and 75th percentiles, within which the medians are indicated by a horizontal line. The bars at the ends of the tails mark the extremes.

CONCLUSION

We formulated a class of nonlinear mixed‐effects statistical models for the estimation of LC100p* from complex studies. Marginal predictions of LC100p*, which support inference over the range of the predictor variables, are easily obtained using only commonly available commercial or open‐source statistical software. We demonstrated six models from that class with the inclusion of two auxiliary predictor variables and two variance components for prediction of effective lampricide formulations from pH. The models outperformed the conventional practice of groupwise estimation of LC100p* followed by secondary modeling of effects of covariates on those estimates by producing more accurate and precise AIC‐weighted averaged predictions of LC50 and LC99.9 in our simulations.

The model class is flexible and is easily extended to accommodate additional features. For example, in other studies, variance might be partitioned across three or more hierarchical levels (e.g., individuals within tanks, among tanks within tests, and among tests). Such additional components of variance require only modification of the specification of the random‐effects formulation. Also, the choice of a function that adequately represents the effect of the covariates is critical to minimization of bias in estimates of LC100p*. That choice should be carefully evaluated for each study.

Fig. 3.

Performance based on precision (root‐mean‐squared error [RMSE]) and accuracy (Bias) of the conventional two‐stage (X̂_2STG) and the marginal mean (X̂_MARG) predictors of (concentrations lethal to 50 and 99.9% of organisms (LC50 and LC99.9, respectively). Performance was estimated from 1,000 Monte Carlo replications of the 38 clusters of the sea lamprey (Petromyzon marinus) trials for which the among‐cluster variance was fixed at σ² = 1.24. Boxes enclose the 25th and 75th percentiles, within which the medians are indicated by a horizontal line. The bars at the ends of the tails mark the extremes. Note that RMSE is scaled logarithmically for visual clarity.

SUPPORTING INFORMATION

Appendix S1. Example SAS^® NLMIXED (Ver 9.1) code for estimation in Equation 11.

Found at DOI: 10.1897/06–630.S1 (23 KB PDF).

Acknowledgements

This work was supported by the Great Lakes Fishery Commission, the U.S. Geological Survey Fisheries and Endangered Aquatic Resources Program in collaboration with the U.S. Fish and Wildlife Service, Hammond Bay Biological Station, and Fisheries and Oceans Canada. We thank Jean V. Adams, Brian R. Gray, and the anonymous reviewers.