A New Stochastic Kriging Method for Modeling Multi-Source Exposure-Response Data in Toxicology Studies.

One of the most fundamental steps in risk assessment is to quantify the exposure-response relationship for the material/chemical of interest. This work develops a new statistical method, referred to as SKQ (stochastic kriging with qualitative factors), to synergistically model exposure-response data, which often arise from multiple sources (e.g., laboratories, animal providers, and shapes of nanomaterials) in toxicology studies. Compared to the existing methods, SKQ has several distinct features. First, SKQ integrates data across multiple sources and allows for the derivation of more accurate information from limited data. Second, SKQ is highly flexible and able to model practically any continuous response surfaces (e.g., dose-time-response surface). Third, SKQ is able to accommodate variance heterogeneity across experimental conditions and to provide valid statistical inference (i.e., quantify uncertainties of the model estimates). Through empirical studies, we have demonstrated SKQ's ability to efficiently model exposure-response surfaces by pooling information across multiple data sources. SKQ fits into the mosaic of efficient decision-making methods for assessing the risk of a tremendously large variety of nanomaterials and helps to alleviate safety concerns regarding the enormous amount of new nanomaterials.

Introduction

Nanomaterials (NM) are finding wide applications in areas such as the energy,[1,2] environment,[3−5] and biomedical engineering sectors.[6,7] The rapid introduction of engineered NM raises imperative concerns on the potential hazard/risk of NM. In comparison to traditional materials and chemicals, it is particularly challenging to fully assess the risks associated with all NM, mainly due to the tremendously large variety of NM. With limited resources for conducting expensive and time-consuming biological experiments, there is an urgent need to develop efficient decision-making methods for NM risk assessment. Primarily motivated by such a need, this paper proposes a new statistical method, which is able to derive more accurate hazard-related information by synergistically modeling multi-source toxicology data. By making more efficient use of given data, the proposed method will help to reduce the experimental cost/time in toxicology studies and contribute to the use of NM in a safe and sustainable manner.

One of the most fundamental steps in assessing the risk of a nanomaterial (or any substance) is to understand and properly characterize its exposure–response relationship.[8,9] A exposure–response relationship describes how the adverse bioactivity effects (the responses) are functionally related to the condition of exposure to a substance.[10] With a well established exposure–response profile, prediction of hazard can be made for a given level of exposure; such a profile also allows for the estimation of exposures at responses of different severities (e.g., the benchmark dose[11,12]), which assists the risk assessor to make judgments to protect a population from increasingly severe effects.

To quantify exposure–response relationships, biological experiments need to be performed under different exposure conditions to observe the corresponding bioactivity responses of animals. Herein, the exposure condition is typically specified through the settings of two quantitative factors: dose level of the substance of interest and time factor involved. Depending on the time scope of the toxicology study, the time factor could be exposure time for long-term studies or post-exposure time for acute studies. On the basis of the experimental data collected, statistical methods are then used to fit exposure–response models quantifying the relationships of interest.

However, it remains a challenge to achieve the exposure–response models of the highest quality from given toxicology data. Existing statistical models[13−15] are not adequate to model typical exposure–response data due to the following three major reasons. (i) First, the exposure–response data for a nanomaterial often arise from multiple sources. More specifically, the data available for modeling typically consists of a number of subsets obtained from different sources such as laboratories,[16] animal providers,[17] and shapes of nanomaterials.[18] The observed response of animals depends not only on the exposure condition but also on these source factors. The subset of data collected from a different source may well reflect a different exposure–response profile. (ii) Second, the exposure–response relationship may well be nonlinear.[19−23] For two-dimensional dose–response curves, a range of nonlinear regression models (e.g., power and logistic models) has been developed[13] to model the sole source data, as opposed to the multi-source data described above. However, the three-dimensional dose–time–response surface is far from being adequately investigated at least partly because of the complex nature of the target surface. (iii) Third, typically the available/affordable toxicology data are not only relatively scarce and highly variable but are also subject to variance heterogeneity.[24] The response variability changes with the experimental setup, which is specified by the exposure condition as well as the data source factors. The variance heterogeneity poses significant difficulty in making statistical inferences for the estimated exposure–response model (i.e., quantifying the uncertainty of the model estimates).

To address these challenges, a new stochastic kriging model is developed in this work, which will be referred to as SKQ (stochastic kriging with qualitative factors). SKQ particularly aims at pooling information from all sources of data to obtain the highest-quality models. It is highly flexible and able to accurately approximate practically any continuous response surfaces,[15,25,26] without requiring a preassumed functional form (e.g., logistic model) as traditional nonlinear regression does.[27] SKQ is able to accommodate variance heterogeneity across exposure conditions as well as different data sources and is able to provide valid statistical inference.

In the literature, the majority of exposure–response modeling methods are restricted to single-source data.[28−30] These methods separately model each source of data and make no effort at all to pool information across different sources. The drawback of such no-pooling methods is that for each data source, a large sample size is required to obtain an adequate exposure–response profile. The mixed effects modeling (MEM) method[14,29,31] represents the closest existing approach for information pooling. (A brief review of MEM is given in Review of Mixed Effects Model, Supporting Information.) However, MEM is parametric regression-based and is subject to restrictive assumptions such as the identification of a specific functional form for regression and the prior-assumed common variance structure for random errors across different data sources. Thus, MEM falls short of addressing challenges (ii) and (iii) stated earlier in this section, whereas SKQ is a more general modeling method.

The remainder of this paper is organized as follows. The Statement of the Research Problem section describes in precise terms the research problem of modeling multi-source exposure–response data in toxicology studies. The new modeling method SKQ is detailed in the Stochastic Kriging with Qualitative Factors (SKQ) section. In the Empirical Studies section, through simulation-based empirical studies, SKQ’s modeling efficiency via information pooling is demonstrated, and the advantages of SKQ over two existing methods are also illustrated.

Statement of the Research Problem

In exposure–response studies for a substance, biological experiments are performed in a range of experimental conditions. An experimental condition is defined by the combination of a number of factors, which can be divided into two categories, quantitative and qualitative factors.

Quantitative factors typically include but are not limited to the toxicant dosage administered to an animal and post-exposure time (or the exposure duration of the toxicant). In this paper, the vector x is used to represent the quantitative factors considered.

Qualitative factors mainly include various source factors such as the laboratory where the experiments were conducted, provider for the experimental animals, and shape of the nanomaterials. The qualitative factors are denoted by the vector z.

The experimental condition is specified in terms of the factor vector w = (xT,zT)T. The random response observed from an animal subject at a factor setting w can be generally written aswhere Y(w) = E[(w)] represents the true expected response, and ε(w) is the random zero-mean error that accounts for the variations across animal subjects.

An example of such exposure–response studies is given later in the subsection Case 1: Modeling Multi-Source Dose-Time-Response Data, where the toxicity of TiO2 nanoparticles (NPs) is investigated. In that case, the vector x = (x1,x2)T includes two quantitative factors: dosage of NPs x1 and post-exposure time x2. There is one qualitative factor z, which has two category levels for the shape of TiO2 NPs: short and long nanobelts.

A setting of the qualitative factors z corresponds to a combination category, say c, and defines a subpopulation or a data source. The total of Q subpopulations specified by the settings of z are denoted as {c; q = 1, 2,...,Q}. In the case of TiO2 NPs, there are Q = 2 subpopulations: c1 for short TiO2 nanobelts and c2 for long TiO2 nanobelts. The population is the union of all the subpopulations and is considered as the TiO2 NPs of both shapes in the aforementioned example. For a given subpopulation c, the bioactivity response obtained from an animal subject is expected to be Y(w|c) = Y(x,c), a continuous function of the quantitative factors x, while subjecting to the cross-subject random error ε(w|c) = ε(x,c).

The biological data collected for the toxicity study of a substance are represented aswhere w denotes the ith design point (factor setting at which experiments are performed) of a total of I distinct design points, (w) is the observed response from the jth replication at w, and n(w) is the number of replications at w.

On the basis of the sample data (eq 2), the objective of statistical modeling is to quantify the dependence of the response upon the quantitative and qualitative factors and to provide valid statistical inference regarding the population being investigated (e.g., TiO2 NPs of both shapes). SKQ aims at achieving this objective with high data efficiency and model generality.

Stochastic Kriging with Qualitative Factors (SKQ)

In this section, the SKQ modeling and inference methods are detailed. As an extension from the standard stochastic kriging (SK), which considers quantitative factors only, SKQ models the variability arising from quantitative as well as qualitative factors. For the readers’ convenience, a review of SK is given in Review of Standard Stochastic Kriging, Supporting Information.

Compared to existing kriging-based methods,[15,32,33] SKQ is the first one that models all of the stochastic elements from the extrinsic variability caused by both quantitative and qualitative factors and intrinsic variability across random replications, as will become clearer in the remainder of this section.

SKQ models the dependence of continuous responses upon the factors w = (xT,zT)T, with x = (x1, x2,..,xd)T ∈ and z = (z1, z2,..., z)T including L qualitative factors. Each qualitative factor z has a number of category levels.

The response at a setting w for the jth replication (animal subject) is modeled by SKQ asThe expectation Y(w) is decomposed into the sum of two parts: f(w)Tβ and M(w). Here, f(w) is a vector of known functions of w, and β is a vector of unknown parameters of compatible dimension. Because it has been widely accepted that f(w)Tβ = β0 (that is, just a constant term) suffices for most applications,[15] this work adopts f(w)Tβ = β0 unless stated otherwise. The term M(w) represents a mean-zero stationary Gaussian process, which intends to capture the extrinsic variability, i.e., the variability due to the factors w.

The randomness of ε(w) is referred to as intrinsic variability. The random noise ε1(w), ε2(w),... at a factor setting w is assumed to have mean zero and be independent and identically distributed (i.i.d.) across replications. The error variance Var[ε(w)] is allowed to be dependent on w.

Given the sample data (eq 2), the sample average of the responses at w across the n(w) replications follows asDenoteas the I × 1 vector of sample average responses at the I distinct design points.

Similarly, the vector of sample average errors is denoted aswith ε̅(w) = n(w)−1 ∑ ε(w),i = 1,2,...,I.

Extrinsic Variance Structure

The modeling of the extrinsic variability is performed following the framework proposed by Qian et al.[32] For the factor settings w = (x,z) and w′ = (x′,z′), the covariance of the stationary Gaussian process M(·) takes the formwith δ2 being the variance of the Gaussian process. The correlation Corr[M(w),M(w′)] is decomposed as the product of two parts: and K(x,x′). To enable the estimation of a SKQ model, specific functional forms need to be assumed for both parts.

In eq 6, K(x,x′) represents the correlation across the quantitative settings for a given combination category of the qualitative factors z and models the variability due to quantitative factors. Hence, K(x,x′) plays the same role in SKQ as in SK, and the discussions regarding K(x,x′) in Review of Standard Stochastic Kriging, Supporting Information, can e inherited here. For specific functional structures of K(x,x′), a range of choices are available in the literature (e.g., Santner et el.,[34] Qian et al.[32]), and one of the most popular correlation functions in practice is the exponential correlation functionIn eq 7, θ = (θ1, θ2,...,θ) is a vector of unknown parameters. It is required that θ > 0 (h = 1,2,...,d), and θ determines the roughness of the response surface for a given combination category of z. The parameter p ∈ (0,2] also needs to be estimated unless p is pre-specified as 2, which corresponds to the widely used quadratic correlation function.[35]

The term in eq 6 is devoted to the correlations across different levels of qualitative factors. As noted in Qian et al.,[32] measures the correlation (similarity) at any two settings w and w′ that differ only on the values of the th qualitative factor. For , a range of functional forms have been proposed in Qian et al.[32] and Zhou et al.[36] Below, two specific correlation functions are given as examples.

• Isotropic (or exchangeable) correlation functions (EC):In eq 8, Φ = {ϕ(; = 1,2,...,L} represents the set of unknown parameters to be estimated, and I[A] is an indicator function that takes 1 if event A is true and 0 otherwise. Clearly, EC assumes that all the category levels of the th qualitative factor are of isotropic nature; that is, a given , is a constant as long as .

• Multiplicative correlation functions (MC):The unknown parameter set Φ includes the following components:c denotes any one of all the possible category levels for the th qualitative factor.

For a given , MC allows the correlation to be dependent on the category levels involved (i.e., z and ).

Given the data collected at I distinct design points, the I × I variance–covariance matrix ∑ is constructed asIn eq 10, R(θ,Φ) denotes the correlation matrix; each component of the matrix represents a correlation, which can be decomposed into two parts as explained above and which is a function of θ and Φ. For an arbitrary setting w0, the I × 1 vector ∑M(w0,·) is defined aswhere v(w0,θ,Φ) denotes the correlation vector with each component being a correlation function dependent on w0, and the unknown parameters θ and Φ.

Intrinsic Variance Structure

The intrinsic variance of the random response at w is denoted as Var[ε(w)], which is dependent on the setting w. Let ∑ε be the I × I variance–covariance matrix of vector ε̅, which is defined in eq 5. Under the i.i.d. assumption for random errors, ∑ε is a I × I diagonal matrix

Estimation and Inference by SKQ (Stochastic Kriging with Qualitative Factors)

The SKQ-based estimation and inference requires the following assumption, which parallels Assumption 1 for SK (see Review of Standard Stochastic Kriging, Supporting Information).

Assumption 1: The random field M is a stationary Gaussian random field; and ε1(w),ε2(w),... are i.i.d. N(0,Var[ε(w)]), independent of ε(w′) for all j and w ≠ w′, and independent of M.

Given a data set and under Assumption 1, Y̅ = (Y̅(w1), Y̅(w2),...,Y̅(w)) as defined in eq 4 follows a multivariate normal distribution with constant mean vector (β0,β0,...,β0)T and variance–covariance matrix

Recall that ∑ and ∑ε are defined in eqs 10 and 12, respectively. Thus, the log-likelihood function of Y̅ in terms of the unknown parameters (β0,δ2,θ,Φ) can be written aswhere 1 is a (I × 1) vector of ones. Following the framework of SK estimation (see Supporting Information), the procedure is developed as follows to obtain the SKQ parameter estimates for (β0,δ2,θ,Φ) that maximize eq 14.

1. Obtain the estimated ∑εwhere

1. Replace ∑ε by and maximize the log-likelihood function (eq 14) with respect to (w.r.t.) (β0,δ2,θ,Φ). Specifically, two steps can be taken to solve the maximum likelihood problem. (i) Given δ2, θ, and Φ, the maximum likelihood estimate (MLE) of β0 is

(ii) Substituting β̂0(δ2,θ,Φ) into eq 14, the problem reduces to maximizingw.r.t. (δ2,θ,Φ), which can be solved by a nonlinear optimization algorithm such as the Matlab fmincon function

1. For an arbitrary setting w0, estimate the expected response Y(w0) bywhere (β̂0,δ̂2,θ̂,Φ̂) are the maximum likelihood estimates obtained from the previous step. Recall that v(w0,θ̂,Φ̂) is defined in eq 11. Following the proof scheme of Ankenman et al.,[15] it can be shown that eq 18 is the best linear unbiased estimator for Y(w0).

The mean squared error (MSE) is obtained aswhere η = 1 – 1T[δ̂2R(θ̂,Φ̂) + ∑̂ε]−1v(x0,θ̂,Φ̂)δ̂2.

The two-sided 100(1 – α)% confidence interval for Y(w0) can be constructed aswhere z1−α/2 is the 100(1 – α/2)th percentile of a standard normal distribution.

Inverse Estimation and Inference

In this section, we consider in particular the SKQ-based modeling/inference for two-dimensional dose–response data and the subsequent derivation of the BMD (benchmark dose) for the substance of interest. The BMD is the dose that corresponds to a specified level of adverse response called the benchmark response (BMR), plays an important role in setting safety standard, and is of particular interest in toxicology studies.

Following the notations adopted earlier, multi-source dose–response data are collected at factor setting w = (x,z), which includes one quantitative factor x representing the dose level and a number of qualitative factors z with all the possible combination categories being {c;q = 1,2···,Q}. The expected dose–response curve is denoted as Y(x,c) for a subpopulation specified by c.

The BMR can be defined as a relative change in the mean response from the control mean or as an absolute level.[13,37] Either definition can be selected based on the knowledge available regarding the substance’s adverse effects, and the BMR defined in one way can be easily converted to that defined in the other. For illustration, we let the BMR be a preselected absolute response in this work, and the BMD for a subpopulation c is written asHere, Y–1 represents the functional dependence of BMD(c) upon BMR, assuming that the inverse mapping exists.[28]

The collection of dose–response curves {Y(x,c);q = 1,2,...,Q} are modeled by SKQ. To perform the inverse calculation as given in eq 21, numerical interpolation needs to be employed based on the fitted SKQ model. In this work, the cubic spline interpolation recommended by Hastie et al.[38] is used to perform the inverse computation for BMD estimation.

Denote as the estimated BMD from the SKQ model for the subpopulation c. The uncertainty of , represented by Var[], directly relates to the safety standard of the substance being investigated. To estimate Var[], numerical methods have to be resorted to again. The bootstrap resampling method developed by Kirk et al.[39] is adapted to quantify the uncertainty of {;q = 1,2,...,Q} based on the SKQ modeling of given data , and the adapted bootstrapping algorithm is described in Figure 1. This algorithm falls into the category of parametric bootstrap resampling method. There are also nonparametric and semi-parametric counterparts,[40−42] and kriging-based bootstrapping inference remains a research topic under investigation.

Figure 1

Bootstrap resampling algorithm for uncertainty quantification of BMD estimates.

Empirical studies suggest that B = 999 usually suffice as a bootstrap sample size for the construction of confidence intervals (CIs). With the B-fitted SKQ models {Ŷ̂*(w);b = 1,2,...,B}, B BMD estimates can be obtained for a specified BMR and a given subpopulation c

On the basis of bootstrapping estimates in eq 22, the nonparametric method suggested by Davison and Hinkley[43] can be easily employed to estimate the 100αth (α ∈ (0,1)) percentile of . In this context, the percentile estimate is referred to as BMDL, which serves as the lower bound of the one-sided 100α% CI for the BMD. For the subpopulation c, the one-sided CI of is written as [,∞). It is worth noting that the percentile method, which is employed here for the estimation of the lower bound BMDL, is one of the various existing methods for estimating the CI boundaries from a given bootstrap sample, and we refer the interested readers to Efron et al.[44] for more information.

Empirical Studies

Empirical case studies were designed and performed to demonstrate SKQ’s advantages to model multi-source exposure–response data over the existing approach, SK (see Review of Standard Stochastic Kriging, Supporting Information) and MEM (see Review of Mixed Effects Model, Supporting Information) method.

Case 1: A multi-source dose–time–response case is developed to show SKQ’s modeling efficiency by pooling information across multiple data sources. The SKQ results are compared to those provided by SK, which models each source of data separately with no information pooling.

Case 2: A multi-source dose–response case is developed to show that compared to MEM (existing information-pooling method) SKQ is a more general method, which is free of the restrictive assumptions stipulated by MEM.

The empirical studies are based on simulation experiments, i.e., sampling through computer experiments whose outputs mimic real experiment data. Simulation, rather than real experiments, is employed for the following reasons. First, in only a simulation-based study, the true expected response surfaces (i.e., the simulation models) are available to evaluate the modeling results. Second, the resulting estimated model is a random outcome, which depends on the random sample data; hence, a modeling method needs to be evaluated in a statistical manner based on the outcomes of applying it on a large number of randomly sampled data sets, which is impossible with real experiments. These advantages of simulation will become clear in the case studies below.

Case 1: Modeling Multi-Source Dose–Time–Response Data

This case is constructed based on the dose–time–response study of TiO2 nanoparticles (NPs) performed by Porter et al.[45] There are two quantitative factors, x = (x1,x2), with x1 ∈ [0,15] μg representing the TiO2 dosage, and x2 ∈ [1,112] days representing the post-exposure time. There is one qualitative factor for the shape of NPs, which is denoted as z. The variable z has two category levels {c1,c2}: c1 denotes short TiO2 nanobelts and c2 long TiO2 nanobelts. Each category corresponds to a different subpopulation/data source. The vector of all the factors is given as w = (x,z). The response of interest is BAL (bronchoalveolar lavage) PMNs measured in the units of 103/mouse.

Simulation Model

The simulation model, which is used to generate simulation data that mimic real experimental data, is described as follows. The true expected responses for the two subpopulations (short and long nanobelts) are represented as {Y(x,c1),Y(x,c2)}, with specific expressions given by Models S25 and S26 in True Expected Exposure–Response Model for Case 1, Supporting Information. Both Models S25 and S26 take the form of a single hidden layer feed-forward neural network and are estimated from real biological data.[46] The true dose–time–response surfaces are plotted in Figure 2.

Figure 2

True exposure–response surfaces for Case 1.

The true variance models used in the simulation are given as

For a subpopulation c (q = 1,2) and at an exposure level x0, a random response y0 is simulated aswhere ϵ is a random error provided by a standard normal random generator.[47]

Case 1 is designed to compare the modeling efficiency of SKQ and SK. MEM has not been applied to this case for the following two reasons. First, MEM relies on a common nonlinear functional form adequate to model the underlying dose–time–response surface for each subpopulation (data source) c, which is very difficult, if not impossible, to identify for the three-dimensional complex surfaces (Figure 2) in this case. Second, the variance structures of eqs 23 and 24 are different across the two categories, which violates the assumption of common variance structure required by MEM. For details regarding the related MEM assumptions, please refer to Review of Mixed Effects Model of the Supporting Information.

Sampling via Simulation

Two types of data sets are obtained via simulation: an estimation data set (EDS) for model estimation/inference and a validation data set (VDS) that is used to evaluate the quality of the model estimated from an EDS.

For Case 1, an EDS includes a total of 32 distinct design points: 16 points for subpopulation c1 (short nanobelts) depicted as stars in Figure 3a and 16 points for subpopulation c2 (long nanobelts) depicted as stars in Figure 3b. At each design point, Model S25 is used to generate eight i.i.d. random responses; that is, eight replications are assigned to each distinct design point.

Figure 3

Design points in the EDS (estimation data set) and check points in the VDS (validation data set).

The VDS includes a dense grid of 16,912 check points, which are depicted as dots in Figure 3. The collection of all the check points is denoted as , with , denotes the collection of dots in Figure 3a, and is the collection of dots in Figure 3b. Further, (q = 1,2) is divided into a number of subsets: , where represents the collection of check points within the subregion specified by the dose range [k – 1,k). The subregions are shown in Figure 3a and b as alternating white and gray rectangles. At each check point, the true expected response Y(·) is available (Figure 2) to evaluate the models fitted from the EDS.

Applying the Modeling Methods

On an EDS generated following the design as given in Figure 3, both SK and SKQ were applied to model the target response surfaces.

SK has no information pooling ability and fits a separate SK model for each subpopulation solely based on the corresponding subset of data. When applying SK, the exponential correlation function (eq 7) is adopted to capture the extrinsic variability for both subpopulations. Two separate SK models are estimated from the short-nanobelt and long-nanobelt data subsets, respectively. At an arbitrary setting w = (x,c), formulas S13 and S14 were used, based on the fitted SK model associated with c, to obtain, respectively, the point estimate and CI of the expected response Y(w).

To apply SKQ, the correlation functions for the extrinsic variability (eq 6) are specified as follows: the exponential correlation function (eq 7) is adopted for K(x,x′), and the EC (isotropic) correlation function (eq 8) is used for . Fitted from an EDS consisting of two sources of data subsets, the resulting SKQ model can be used for estimation and inference. At an arbitrary factor setting w = (x,c), formula (eq 18) provides the point estimate for Y(w), and eq 20 gives the CI for Y(w).

Pooling Effects of SKQ

SKQ’s strength lies in its ability to pool information across multiple subpopulations/data sources. The estimation quality w.r.t. a subpopulation can be substantially improved because SKQ allows for the borrowing of information (or data) from all the other subpopulations. To demonstrate SKQ’s estimation efficiency, herein we compare the SK and SKQ results in terms of the quality of their point estimates for Y(·).

The estimated root mean squared error (ERMSE) defined as follows is used to evaluate the goodness of point estimate Ŷ(·). Recall that is the collection of the checkpoints for subpopulation c included in the rectangle specified by dose range [k – 1,k), as shown in Figure 3. We definewhere represents the total number of check points in the set . Clearly, eq 26 measures the average deviation of Ŷ(·) from the true value Y(·) at the check points in .

Applying a modeling method (SK or SKQ) on one EDS, denoted as EDS( (r = 1,2,...,R), is considered as one macro-replication and leads to a set of performance statisticsFor empirical evaluation, a total of R = 1000 independent EDS has been generated by simulation experiments, and the average ERMSE across the macro-replications are denoted asand are summarized in Table 1. The statistic measures the expected overall deviation between Y(·) and Ŷ(·) within the subset . As shown in Table 1, in each subset of check points , the point estimates given by SKQ are much more accurate than those provided by SK. For 22 out of the 30 check point subsets in Table 1, the average ERMSE given by SKQ is less than half of that given by SK. Clearly, by synergistically modeling multi-source data, SKQ leads to statistical models of substantially improved quality.

Table 1

Comparison of Estimation Results from SK and SKQ for Case 1

subset of check points	SK	SKQ	subset of check points	SK	SKQ
	1.6110	1.5876		5.1988	2.9626
	2.0351	1.4277		4.1851	2.5193
	2.4907	1.2814		3.4493	2.1025
	2.9825	1.1738		2.8623	1.6950
	3.4742	1.1290		2.4246	1.3363
	3.9355	1.1518		2.1431	1.0676
	4.3396	1.2258		1.9998	0.9487
	4.6247	1.3294		1.9694	1.0057
	4.7185	1.4512		2.0524	1.1582
	4.5792	1.5929		2.2103	1.3311
	4.2379	1.7578		2.3942	1.4815
	3.8478	1.9497		2.5792	1.5902
	3.7447	2.1753		2.7650	1.6665
	4.3577	2.4427		2.9780	1.7489
	5.8143	2.7547		3.2707	1.9007

Case 2: Modeling Multi-Source Dose–Response Data

This case is derived from the dose–response study of TiO2 nanoparticles (NPs) after 3 days of exposure performed by Porter et al.[45] There is one quantitative factor x representing the dose level, and x ∈ [0,20] μg. The one qualitative factor is denoted as z with three categories {c1,c2,c3}. Each category is used to represent a different batch of animals. The vector of factors is written as w = (x,z). The response of interest is BAL (bronchoalveolar lavage) PMNs measured in the units of 103/mouse.

Simulation data that mimic the real experiment data are generated using the following true dose–response modeland variance model

For a certain category c and at a dose level x0, a random response y0 is simulated aswhere ϵ is a random error provided by a standard normal random generator.[47]

This case is used to compare SKQ and MEM and is designed in such a way that the two basic assumptions required by MEM are met: (i) A nonlinear functional form can be easily identified and employed to model the target dose–response curves. (ii) There is a common variance structure (eq 28) across different data sources. The third assumption made by MEM is the multivariate normality of the model coefficient vector (see Review of Mixed Effects Model, Supporting Information). According to eq 27, the three true coefficient vectors in this case are (20,11)T, (19.5,10.5)T, and (26,13.2)T; on the basis of which it is hardly possible to judge whether the normality assumption holds or not. This is quite typical of multi-source data.

Simulation-Based Sampling

As in Case 1, both EDS and VDS are generated in this study for model estimation and evaluation, respectively. To generate an EDS, the design as shown in Table 2 is used; eight replications are carried out at 5 evenly spaced dose levels for each of the three categories. An example EDS is given in Table S1 of An Estimation Data Set for Case 2 of the Supporting Information.

Table 2

Design Points in EDS (estimation data set) for Case 2

x: dose	0	5	10	15	20	0	5	10	15	20	0	5	10	15	20
cq: subpopulation	c1	c1	c1	c1	c1	c2	c2	c2	c2	c2	c3	c3	c3	c3	c3

The two alternative information pooling methods, MEM and SKQ, were applied on the EDS given in Table S1 of the Supporting Information.

To perform MEM, the common nonlinear functional form of the dose–response curve for each of the three subpopulations is assumed to bewith α = (α,α) being the unknown parameters for subpopulation c (q = 1,2,3). The variance model is assumed to followwhere σ and γ are unknown parameters common to different subpopulations. Note that the assumed forms (eqs 30 and 31) are the ones that the true models (eqs 27 and 28) fall into, respectively; the adoption of functional forms (eqs 30 and 31) is meant to allow MEM to achieve its ideal estimation results.

With the assumed forms (eqs 30 and 31), MEM is performed on the EDS in Table S1 of the Supporting Information. The fitted dose–response models areThe fitted variance model isWith the fitted MEM, the expected response can be estimated at any dose level, and the BMD estimate can be obtained for a given BMR. The confidence intervals can also be constructed for these quantities of interest (see Review of Mixed Effects Model, Supporting Information).

When applying SKQ to the same EDS (Table S1, Supporting Information), the correlation (eq 6) is constructed as follows: The exponential correlation function (eq 7) is used to model the correlations between quantitative variables, and the MC correlation function (eq 9) is used to model the correlations across different levels of qualitative factors. Normalization of the original data (Table S1, Supporting Information) was also performed so that both the quantitative factors and responses range over [0,1]. Applying the maximum likelihood estimation procedure (see Estimation and Inference by SKQ) on the normalized data leads to the fitted parameters for the SKQ model as displayed in Table 3.

Table 3

SKQ Parameters Estimated from Normalized Dose–Response Data for Case 2

β̂0	δ̂2	θ̂1	p̂	φ̂1	φ̂2	φ̂3
0.5426	0.1272	1.3588	1.9168	0.01	0.0144	0.022

With the fitted SKQ model, the expected response can be estimated at any dose level (see Estimation and Inference by SKQ), and the inverse BMD can be obtain numerically (see Inverse Estimation and Inference); the confidence intervals for these quantities can also be obtained accordingly. Note that when utilizing the SKQ specified by Table 3 for the estimation/inference, a simple conversion calculation is needed to ensure that the estimates are given on the original scale because those parameters are obtained from the normalized data.

Comparison of the Two Modeling Methods

The modeling results of MEM and SKQ are compared in terms of their estimation/inference abilities for (I) the expected responses as well as (II) the BMD values.

(I) Estimation/Inference of the Expected Response Ŷ(·)

From the one EDS given in Table S1 of the Supporting Information, both MEM and SKQ were applied as described earlier, and the estimation results for the two methods are displayed in Figure 4. The circles denote the EDS and are plotted in both Figure 4a and b. The dashed curves represent the estimated expected responses, and the solid curves are the lower and upper 95% CI curves for the true expected responses. As shown in Figure 4, over the dose range, the widths of the CIs (that is, the vertical distances between the lower and upper CI curves) provided by SKQ are generally narrower than those given by MEM; this pattern, which is graphically illustrated for the one EDS plotted as circles in Figure 4, holds consistently for all of the R = 1000 macro-replications carried out in this study.

Figure 4

Comparison of the dose–response fitting results from MEM and SKQ.

As explained in Case 1, one macro-replication refers to the process of applying a method on one randomly generated EDS. Using a modeling method (MEM or SKQ), R = 1000 macro-replications lead to 1000 CIs of the true expected response Ŷ(·) for any check point specified in terms of (x,c). Hence, the coverage probability of the CIs can be estimated as the percentage of the 1000 CIs that include the true expectation Ŷ(·). In our simulation study, Ŷ(·) is available from the simulation model (eq 27) for the purpose of evaluating the CIs. Ideally, among these 1000 CIs, the percentage of the CIs that actually contain Ŷ(·) should be very close to 95%, the nominal coverage level.

Table 4 presents the coverage probabilities of the 95% CIs given by MEM and SKQ, respectively, based on each method’s 1000 macro-replications. The first two rows of Table 4 specify a number of check points. The estimated coverage probabilities of the MEM CIs are given in the row marked as “MEM’’, and are all 1.000 at the check points, which is much higher than the nominal 95%. The estimated coverage probabilities resulting from SKQ are given in the row labeled as “SKQ’’ and are much closer to the nominal percentage 95%.

Table 4

Comparison of MEM and SKQ in Terms of CI Coverage Probabilities for the Expected Response Y(·)

	subpopulation C1				subpopulation C2				subpopulation C3
x:dose	2.5	7.5	12.5	17.5	2.5	7.5	12.5	17.5	2.5	7.5	12.5	17.5
MEM	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000
SKQ	0.957	0.966	0.974	0.933	0.952	0.955	0.970	0.926	0.955	0.970	0.964	0.968

Therefore, as shown in Figure 4 and Table 4, SKQ is able to provide tighter CIs (i.e., CIs that are narrower and with more on-target coverage probabilities) for the true expected responses, while MEM overshoots the nominal coverage percentage by providing wider CIs.

(II) Estimation/Inference of the BMD

Herein, the two methods, MEM and SKQ, are compared in terms of their inverse estimates for the BMD associated with a pre-specified BMR. For demonstration, the BMR is set as 42 × 103/mouse in this case. (Recall that the response of interest is BAL PMNs measured in the units of 103/mouse).

As already explained, R = 1000 macro-replications were performed using each of the two methods based on the 1000 data sets {EDS(;r = 1,2,...,R}. From the rth (r = 1,2,...,R) macro-replication, a one-sided 95% CI was constructed for BMD(c), q = 1,2,3, following the bootstrapping resampling method (Inverse Estimation and Inference); the lower bound of the one-sided CI is called BMDL. The BMDLs estimated from the 1000 MEM macro-replications are denoted asThe BMDLs obtained from the 1000 SKQ macro-replications are represented asThe true values for BMD(c), q = 1,2,3, can be easily obtained for BMR = 42 × 103 based on the simulation model (eq 27) and are represented by the horizontal lines in Figure 5a–c corresponding to the three subpopulations. Figure 5a is devoted to the subpopulation c1, and the two box plots are generated based on {r = 1,2,...,1000} and {r = 1,2,...,1000}. Figure 5b and c are plotted for the other two subpopulations in the same manner. Clearly, the BMDLs given by SKQ are much closer to the true BMD than those provided by MEM.

Figure 5

Box plots for the BMDLs resulting from the two modeling methods.

Each BMDL estimate in eqs 34 and 35 corresponds to a one-sided 95% CI: [BMDL,∞]. Table 5 compares the coverage probabilities of the CIs obtained from MEM and SKQ. It is shown that the coverage probabilities of SKQ are close to the nominal level 95%, whereas MEM’s estimated probabilities are all 1.000. Therefore, as shown in Figure 5 and Table 5, SKQ is able to give more informative CIs of the BMD compared to MEM.

Table 5

Comparison of MEM and SKQ in Terms of CI Coverage Probabilities for BMD

	subpopulation C1	subpopulation C2	subpopulation C3
pre-specified BMR	42	42	42
MEM	1.000	1.000	1.000
SKQ	0.9349	0.933	0.959

Summary

A new semi-parametric statistical model, SKQ (stochastic kriging with qualitative factors), has been developed in this work. SKQ is the first kriging-based model able to take into account the three types of variability steming from quantitative factors, qualitative factors, and uncountable sources (random errors). Compared to the parametric MEM (mixed effects modeling) method, the closest regression-based counterpart, SKQ represents a more general modeling approach and is free of the various restrictive assumptions stipulated by MEM.

Through the empirical simulation studies, the modeling efficiency of SKQ is demonstrated over the existing methods, SK (standard stochastic kriging) and MEM. SKQ is able to pool information across multiple data sources, accommodate general data features, and provide more informative estimation/inference from given data. For clarity and succinctness of the presentation, the two cases for this paper were designed to involve a relatively small number of data sources (subpopulations). It is worth noting that when applying SKQ to data of a large, as opposed to a small, number of sources, there is hardly any additional theoretical or implementation hurdles. In our empirical experience, SKQ’s modeling efficiency (information pooling effects) is more pronounced with more sources of data. All the SKQ-related algorithms are currently implemented in MATLAB programs, which are available upon request from the authors. To facilitate SKQ modeling and inference by biologists, a VBA (Visual Basic for Application) Excel macro tool is under development. The VBA tool intends to allow users who possess the basic ability to use Excel to perform SKQ by interacting with Excel spreadsheets and menus.

As a final note, the focus of this work is on efficient modeling of given data. However, SKQ’s ability to provide more informative statistical inference from limited data certainly opens opportunities for efficient design of experiments (DOE). Given some subsets of data already collected and likely from different sources, how should we efficiently design the next stage biological experiments so that all the integrated data is most informative? DOE is directly associated with model estimation and inference, and in our ongoing research, a SKQ-based DOE method is under development to achieve experimental efficiency.