Estimating equations for biomarker based exposure estimation under nonsteadystate conditions
 Scott M Bartell^{1, 2}Email author and
 Wesley O Johnson^{2}
DOI: 10.1186/1476069X1057
© Bartell and Johnson; licensee BioMed Central Ltd. 2011
Received: 8 October 2010
Accepted: 13 June 2011
Published: 13 June 2011
Abstract
Unrealistic steadystate assumptions are often used to estimate toxicant exposure rates from biomarkers. A biomarker may instead be modeled as a weighted sum of historical timevarying exposures. Estimating equations are derived for a zeroinflated gamma distribution for daily exposures with a known exposure frequency. Simulation studies suggest that the estimating equations can provide accurate estimates of exposure magnitude at any reasonable sample size, and reasonable estimates of the exposure variance at larger sample sizes.
Background
Health risks assessment, dietary research, environmental epidemiology, and other endeavors that depend on quantitative exposure estimation are increasingly making use of exposure biomarkers instead of, or in addition to, traditional contactbased exposure estimates [1–4]. Estimation of exposure rates from biological measurements is a particularly challenging problem given the complex relationship between ingested, inhaled, or dermally absorbed chemical exposures and the resulting tissue concentrations over time. Indeed, because many different exposure patterns can lead to the same blood or urine concentration at a given point in time, typical approaches to biomarker based exposure estimation must rely on simplifying assumptions regarding the exposure patterns. In cases where exposures have ceased, such as postshift or postretirement studies of occupational exposures, exposure rates are often reasonably assumed to be zero. In other settings, investigators often rely on an assumption that exposure rates are constant over time for each individual.
Unfortunately, virtually all environmental exposures vary in magnitude over time, thereby violating the steadystate model. For example, ingestion occurs intermittently and only during waking hours. These violations can cause substantial errors in biomarker based exposure estimates that rely on a steadystate assumption [5]. The degree of error introduced by the steadystate model is often substantial depending on the elimination rate of the chemical, the frequency of contact, and the variability in exposure over time, even under the highly optimistic assumption that every individual's exact biokinetic parameters are known [5–7]. These results suggest that a substantial portion of observed population variability in mercury biomarker concentrations may result from nonsteadystate exposure conditions, rather than being entirely attributable to true differences in individual mercury exposure rates. Exposure measurement error is known to cause bias in epidemiologic dose response modeling, though posthoc methods of adjustment have been proposed [8].
Formal statistical methods for biomarker based exposure estimation that do not rely on steadystate assumptions are needed. Standard Monte Carlo simulation methods have been suggested but are inadequate for inverse estimation problems due to unknown but nonzero correlations [9].
We present a new statistical method for estimating individual exposures to mercury based on individual hair or blood mercury biomarkers and individual exposure frequencies, for a group of people with the same probability distribution of daily exposure magnitudes. Although this method was developed using mercury as a case study, it may be applicable to other toxicants. This method avoids steadystate assumptions and incorporates information from biokinetic models. The new approach utilizes discretetime approximations to continuoustime biokinetic models and statistical methods based on the theory of estimating equations [10, 11].
Methods
Simplified Biokinetic Model
The relationship between chemical exposure at the visible exterior boundary of a person and related biomarker measurements is governed by fairly complicated timedependent processes including absorption across the skin, gastrointestinal tract, or lung epithelium; distribution throughout the body via blood circulation, filtration, metabolism, and/or sequestering by liver and kidneys; and excretion via skin, nails, hair, urine, and feces. The entire system of these processes is referred to as "pharmacokinetics," "toxicokinetics," or "biokinetics."
Mathematically, biokinetic models are typically composed of continuoustime systems of differential equations, with each differential equation representing the rate of change in concentration or mass of a chemical in a particular tissue or organ as the chemical is exchanged with blood, metabolized, or excreted. These biokinetic models are generally not invertible without additional constraints and ad hoc methods, due to the dimensional reduction from continuoustime exposure patterns to biomarker measurements at specific time points. In other words, many different exposure patterns can lead to the same biomarker concentration, so it is generally not possible to determine an exact exposure pattern using only biomarker measurements. Instead, biomarker based risk assessments typically rely on the simplifying but unrealistic steadystate assumption, multiplying each individual's biomarker concentration by a steadystate ratio in order to estimate constant exposure rates.
The classic single compartment biokinetic model can be expressed as a differential equation: , where y(t) is the biomarker concentration at time t, f is the fraction of ingested mercury present in the blood after absorption across the gastrointestinal tract and equilibration throughout the body, v is the volume of blood, I(t) is the mercury exposure rate at time t, and k is the excretion rate coefficient for mercury. If I(t) is constant then y(t) eventually reaches "steadystate," i.e., and . Thus, the steadystate ratio of I to y is . An estimate of this steadystate ratio is typically multiplied by each individual's measured biomarker concentration in order to estimate the corresponding exposure rate.
where y_{ it }is the blood mercury concentration in individual i on day t, I_{ ij }is the mercury intake for individual i on day j, W_{ ijt }is the "weight" or influence of the day j intake on the day t biomarker measurement in individual i, and ε_{ it }is a statistical error term with expectation 0 and variance σ_{ ε }^{2}. W_{ ijt }expresses the short term absorption and dilution of the mercury into the blood, as well as the elimination of mercury from the body over time. As noted by Sherlock et al. [12] in their least squares estimation of biokinetic parameters from a controlled dosing study, Equation 1 provides a close approximation to a continuoustime onecompartment biokinetic mercury model with first order elimination, provided that , where f_{ i }is the fraction of ingested mercury present in the blood of individual i after absorption across the gastrointestinal tract and equilibration throughout the body, v_{ i }is the volume of blood in individual i, and k_{ i }is the excretion rate coefficient for mercury. For chronic exposures t should be chosen to reflect at least ten halflives, but needn't include the entire individual history as the earliest exposures will have negligible contributions to the measured biomarker concentration.
where h_{ i }is the equilibrium ratio for hair to blood in individual i [5, 13].
Probability Model for Daily Exposures
where ω_{ i }is the exposure frequency (with units of day^{1}), 1_{{S} }is an indicator function that takes the value 1 when statement S is true and 0 otherwise, and a and λ are parameters describing the gamma distribution. Unlike the lognormal distribution, the gamma distribution can take on a heavily skewed shape or a nearly symmetric shape depending on the values of the two parameters. Here we assume that exposures are independent across days and across individuals. This independence assumption may not be reasonable for individuals who share meals, or for those who obtain multiple meals from the same source item. For example, a person consuming many tuna steaks all cut from the same individual fish should have highly correlated mercury exposures over time.
One important attribute of the zeroinflated gamma distribution shown in Equation 3 is that its expectation and variance are easily computed: and . For our proposed estimation method, it is useful to reparameterize the zeroinflated gamma distribution using and . α and β are the log variance and log mean of the conditional exposure distribution, for days with nonzero exposures. The expectation and variance of the unconditional exposure distribution can then be written as and . We will estimate α and β rather than a and λ. This parameterization has two important advantages: 1.) α and β are unrestricted on the real number line and 2.) E(I_{ ij }) has only one unknown parameter when individual exposure frequencies can be measured, e.g. by food frequency questionnaires.
In this model all methylmercury exposures within a day are grouped together, so the exposure frequency cannot exceed 1 per day. Although the daily grouping of exposures represented by Equations 13 does not capture the full complexity of the exposure profile, the approach is much more realistic than the assumption of constant mercury exposure, is amenable to formal statistical treatment, and can easily be extended to include fixed covariate effects.
We have chosen to group exposures by day, but our model is easily adjusted for grouping into smaller (or larger) time intervals, provided that the biokinetic model weights are selected appropriately. For the best approximation, interval lengths should be small relative to the biological half life of the toxicant being modeled. For example, many solvents are quickly excreted from the body; one day exposure aggregates for these compounds would be too crude to compare with biomarker concentrations, but one to sixty minute intervals might prove reasonable.
Results
In the case of exposure estimation using biomarker measurements alone, the models formulated by combining the above equations pose a challenge in that the likelihood equations are difficult to obtain due to the convolution of many mixture distributions containing both discrete and continuous components. We first explored normal approximations to the summation in Equation 1, but simulation studies indicate that normality only holds when the exposure frequency is high and the variance in daily exposure magnitudes is low, making the normal approximation and classical statistical methods unsuitable for most realistic exposure settings [13]. Instead, we propose estimating equations based on the quasilikelihood [11]. The estimating equations rely entirely on the expectation and variance of the biomarker measurements in terms of the unknown exposure parameters, bypassing the need for an explicit likelihood equation or even specification of exact probability distributions.
Algorithm
Estimating equations, particularly in the form of generalized estimating equations [14], have become popular in situations where it is difficult to model complex data, such as correlated data that do not arise from a multivariate normal distribution [11, 15]. Although estimating equations do not appear to have been previously applied to nonsteadystate biomarker based exposure estimation, the method is quite flexible and appears to be reliable in this setting.
These methods make use of a concept called the "quasiscore function" [10, 16]. Consider an n × 1 response vector Y with expectation vector EY and covariance matrix V. Let EY be a function of an unknown pparameter vector β, and D be the n × p matrix . The quasiscore function is the p × 1 vector.
starting with , an initial guess for β. , , and E_{1}Y are the lth iterate estimates of D, V, and EY, respectively, and are all obtained by evaluation at . The algorithm consists of repeated application of Equation 4, incrementing l by 1 each time, until and only differ by a prespecified negligible amount, at which point has usually converged to the root of the estimating equations. is an estimate of the covariance of , and is easily obtained directly from the algorithm. Although we do not show it here, Equation 4 can usually be written in a more computationally efficient form involving summations of block diagonal elements [15].
In a simple biomarker application, Y might represent statistically independent blood mercury measurements y_{ it }for i = 1, 2, ..., n individuals, with one measurement per person. Assume that Equations 13 apply and that each individual's biokinetic parameters k_{ i }, f_{ i }, and v_{ i }and exposure frequency ω_{ i }are known. Using Equations 13, basic mathematical properties of expectations and variances of weighted sums [17], and summation rules for finite geometric series, the vector EY consists of the n elements , i = 1, 2, ..., n, and V is a diagonal matrix with diagonal entries . [13]. In this case, p = 1 because there is only one unknown parameter β in the mean vector EY. Conveniently, here D = EY.
In most cases, including our model for biomarkers, V depends on additional unknown parameters other than β. These additional parameters are denoted by αa scalar in our model, as there is only one unknown variance parameter not contained in the mean vector. There are several different strategies for estimating both α and β, but the most reliable appears to be the use of alternating estimating equations, whereby a second estimating equation is written for α, and the algorithm proceeds with alternating iterative estimation of α and β [13, 15]. The estimating equation for α can be written as , where is the upper diagonal of the estimated covariance matrix in vector form as , s is an "empirical covariance vector" with (n^{2}+n)/2 elements s_{ ij }= (Y_{ i } EY_{ i })(Y_{ j } EY_{ j }) corresponding to the elements of , is the (n^{2}+n)/2 length vector of estimates for , and is the (n^{2}+n)/2 x (n^{2}+n)/2 estimated covariance matrix for the vector s.
When the observations are independent (a reasonable assumption for one measurement per person), V is a diagonal matrix and many of the elements of are 0. In this special case, the estimating equations for α can be simplified using with a corresponding n length vector for s, an n length vector for , and an n x n matrix for . In either case, an iterative equation analogous to Equation 4 can be derived from the estimating equations for α. Because there is only one measurement per person, this method relies on the betweensubject variability in biomarker measurements for estimation of α. This approach is reasonable when subjects have similar exposure sources (e.g., similar types of fish in the diet).
The elements of the vector D * are easily obtained as . An expression for V * is more difficult to obtain without simplifying assumptions. We employ the "independence working matrices" assumption of Prentice and Zhao [15], approximating the elements of V * based on a simplifying assumption of independence and normality among the elements of Y. For one measurement per person, this assumption results in a diagonal matrix for V *, with elements (2V_{11}, 2V_{22}, ..., 2V_{nn}). In the case of multiple measurements per person, V * is block diagonal with covariance terms of V_{ ii }V_{ jj }for measurements at different time points in the same individual. Though these can be crude approximations for V *, resulting estimates of α and β remain theoretically valid and appear to be reliable at reasonable sample sizes.
It is impossible to estimate both α and from the data alone with only one biomarker measurement per person; in this setting an external estimate of one of the two parameters is needed. In other words, one needs to provide an estimate of either 1.) the variance in exposure magnitude over time, or 2.) the variance of any random errors resulting from sources other than exposure variability (e.g., biomarker measurement error and natural biokinetic variation across individuals. Collecting multiple biomarker measurements per individual may reduce or eliminate the need for specifying external parameter estimates for α or .
Testing and Implementation
In order to examine the performance of the estimating equations under simple conditions, we conducted several simulation studies for the simple setting described above, with only one blood measurement per individual. We generated 10,000 data sets at each of 9 different exposure frequencies and 10 different sample sizes ranging from 2 to 1024. For simplicity we assigned all individuals the same biokinetic parameters throughout the simulations, using the values f_{ i }= 0.0475, v_{ i }= 5 L, and k_{ i }= 0.014 d^{1} [1, 18] for all i, and assumed that σ_{ ε }^{2} was known and relatively small. For ease of interpretation we present results in terms of the mean exposure magnitude, μ = e^{ β }, and the variance in exposure magnitudes, . All simulations were performed using μ = 10 μg d^{1}, σ_{ g }^{2} = 5 μg^{2} d^{2}, σ_{ ε }^{2} = 0.03^{2} μg^{2} d^{2}, and t = 1000 d. The algorithm only failed to converge for a few simulated data sets with both low exposure frequency and a sample size less than 10; for the worst case with n = 2 and an exposure frequency of 0.1 the algorithm converged for about 99.4% of the simulated data sets.
Bias in estimation of the mean exposure magnitude, μ
Exposure Frequency  

Sample Size  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9 
2  0.033  0.024  0.021  0.008  0.007  0.003  0.002  0.005  0 
4  0.002  0.007  0.017  0.002  0.006  0.007  0.002  0.004  0.005 
8  0.01  0.004  0.004  0.013  0.001  0.001  0.001  0.003  0.001 
16  0.008  0.001  0.003  0.007  0  0  0.001  0.001  0.001 
32  0.001  0.004  0.003  0.007  0.001  0.003  0.001  0.002  0.001 
64  0.007  0.004  0.003  0.005  0.003  0  0  0.001  0.001 
128  0.005  0.004  0.002  0.007  0.002  0.001  0  0.001  0.001 
256  0.006  0.003  0.001  0.005  0.001  0.001  0.001  0.001  0.001 
512  0.004  0.001  0.002  0.007  0.002  0.001  0.001  0.001  0.001 
1024  0.004  0.002  0.002  0.007  0.002  0  0  0.001  0 
Bias in estimation of the variance in exposure magnitudes, σ_{ g }^{2}
Exposure Frequency  

Sample Size  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9 
2  5.34  5.9  4.8  3.36  2.89  1.27  0.3  0.82  1.83 
4  10.22  8.55  6.89  6.05  4.74  3.11  1.96  0.26  0.96 
8  9.22  8.49  6.89  5.48  4.56  2.91  1.75  0.41  0.56 
16  7.6  6.31  5.29  4.45  3.07  2.38  1.07  0.09  0.4 
32  5.36  4.56  3.44  2.56  2.31  1.43  0.54  0.02  0.35 
64  3.38  2.88  2.21  1.69  1.39  0.92  0.17  0.12  0.17 
128  2.09  1.58  1.32  0.77  0.68  0.36  0.03  0.19  0.06 
256  0.95  0.72  0.61  0.44  0.31  0.18  0.1  0.09  0 
512  0.43  0.31  0.2  0.14  0.22  0.21  0.13  0.09  0.02 
1024  0.05  0.06  0.03  0.01  0.19  0.16  0.08  0.04  0.04 
Coverage of nominal 95% confidence intervals for the mean exposure magnitude, μ
Exposure Frequency  

Sample Size  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9 
2  94.5%  94.7%  94.8%  94.8%  94.3%  94.0%  94.1%  93.7%  91.0% 
4  94.9%  95.0%  94.8%  95.0%  94.8%  94.6%  94.7%  93.8%  91.9% 
8  95.3%  95.5%  94.9%  95.2%  94.9%  94.4%  94.6%  94.4%  93.2% 
16  95.3%  95.3%  95.2%  94.9%  95.2%  94.9%  94.8%  94.0%  93.4% 
32  95.1%  95.5%  95.3%  95.1%  95.1%  94.7%  94.8%  94.4%  94.2% 
64  95.1%  95.2%  95.7%  94.8%  94.9%  94.5%  94.5%  94.5%  94.3% 
128  95.2%  94.9%  95.1%  95.1%  95.2%  94.6%  94.8%  94.4%  94.9% 
256  95.2%  95.1%  95.2%  94.9%  95.1%  95.2%  94.9%  95.0%  94.9% 
512  95.5%  95.3%  95.2%  94.5%  94.9%  94.8%  94.8%  95.0%  95.4% 
1024  95.1%  94.9%  95.0%  94.3%  95.2%  95.2%  94.9%  94.9%  94.7% 
Discussion
The simulation studies indicate that the estimating equations are quite reliable for estimation of mean exposure magnitudes even at fairly low sample sizes, but that the variance in exposure magnitudes may be difficult to estimate at low sample sizes when the exposure frequency is low. Confidence intervals are also readily obtained and reasonable accurate for mean exposure magnitudes, but may be overly conservative for the variance in exposure magnitudes.
It is worth noting that the steadystate method assumes a constant exposure rate, implying that the exposure frequency is 100% and that the exposure variance is 0. Observed betweenindividual variability in biomarker measurements is therefore assumed to reflect true individual differences in exposure rates, rather than daytoday variability in exposures. The steadystate method has previously been shown to be imprecise in estimating individual exposure rates when daytoday variability exists [5]. If the correct value is used for the steadystate ratio b, the steadystate method estimates μ with a bias of and a precision of when applied to an individual with nonsteadystate zeroinflated exposures [13]. In our simulation setting, the bias of the steadystate estimate for μ is therefore approximately 10(ω_{ i } 1) and the variance of the steadystate estimate (for each individual) is approximately . Thus, the standard errors for the steadystate estimate range from 0.26 to 0.44 in our simulation setting. Although the estimating equations appear to produce unbiased estimates for μ in nearly all of our simulations, clearly outperforming the steadystate estimate at all nine exposure frequencies in this setting, it is worth noting that a simple modification of the steadystate estimator from to produces nearly unbiased estimates for large t, and averaging those estimates for groups of individuals with similar exposure sources would improve the precision of the estimate of μ. With these modifications, the steadystate method might be a reasonable approach for estimating groupaveraged exposure rates when the exposure duration is long, provided that exposure variability is only a nuisance instead of a target for estimation.
Unlike steadystate methods, the estimating equations provide estimates and standard errors for both exposure magnitude parameters when σ_{ ε }^{2} is negligible or can be estimated from external data. Statistical theory and our simulations suggest that the estimating equations estimates of β have approximately normal distributions at sample sizes above 20 or 30. At high exposure frequencies, even fewer biomarker samples may be necessary in order for the estimates to be normal. Unfortunately, estimates of α generated in the simulation studies exhibited fairly strong departures from normality even with hundreds of samples, suggesting that the usual asymptotic normal confidence intervals for α might not be appropriate for typical biomarker studies. If confidence intervals for α are desired, jackknife or bootstrap procedures might provide more accurate results. However, such estimates still depend on the correct specification of σ_{ ε }^{2}. If variability in daily exposure magnitudes is an important target for inference, we recommend that multiple biomarker measurements be obtained for each individual. Carefully structured repeated biomarker measurements and duplicate samples may provide a means to simultaneously estimate both α and σ_{ ε }^{2}.
Further study is needed to determine whether this approach or alternatives such as populationaveraged generalized estimating equations are more reliable in this setting.
Incorporation of interindividual variability in the biokinetic parameters is another goal for extension of these methods. Although additional variance parameters would be difficult to estimate with only one measurement per person, biokinetic parameters that vary across individuals but are stable over time might be estimable from repeated biomarker measurements. In principle these extensions might be accomplished using estimating equations, but Bayesian approaches become more attractive with increasing model complexity [13].
Conclusions
Direct exposure measurements such as those obtained by duplicate diet studies are often prohibitively expensive for chronic exposure situations such as mercury exposure via ordinary seafood consumption. Indirect estimation using exposure biomarkers is an informative and less expensive approach, but such an exercise should be recognized as a statistical problem whereby the unknown exposure parameters are estimated based on a theoretical model relating the unknown exposures to the observed biomarker measurements.
Our proposed estimating equation approach to biomarker based exposure assessment represents a compromise between the steadystate model, which is overly simplistic but still widely used because of its practicality, and fully detailed biokinetic models that are somewhat impractical for use in formal statistical estimation with ongoing exposures. Our methods have some clear advantages for mercury exposure estimation compared to the steadystate model, due to the more realistic model and the ability to do hypothesis testing and statistical inference. We also believe it may be a valid approach for other chemicals exhibiting firstorder biokinetics, provided that the discretetime unit length is selected to be short relative to the biological halflife. However, its current implementation is based on somewhat restrictive assumptions, including that biokinetic parameters are known and constant across individuals, that measurement error is negligible or can be estimated externally, that individuals can be grouped according to similar exposure distributions, that exposures are independent across days and individuals, and that individual exposure frequencies can be accurately measured.
Future work should assess the performance of both steadystate and nonsteadystate methods when these assumptions are violated, as well as extending these methods towards less restrictive assumptions.
List of Abbreviations
 EPA:

Environmental Protection Agency
Declarations
Acknowledgements
This research was supported by a graduate student fellowship for SB from the US Environmental Protection Agency's Science to Achieve Results (STAR) program. This manuscript has not been subjected to EPA review, and no official endorsement should be inferred. The authors would also like to thank Dr. Laurel Beckett for suggesting generalized estimating equations as a potential approach for this problem.
Authors’ Affiliations
References
 Environmental Protection Agency [EPA]: Integrated Risk Information System: Methylmercury. [http://www.epa.gov/iris/subst/0073.htm]
 Ryan L, Huang W, Thurston SW, Kelsey KT, Wiencke JK, Christiani DC: On the use of biomarkers for environmental health research. Stat Methods Med Res. 2004, 13: 207225.
 Lu C, Barr D, Pearson M, Bartell S, Bravo R: A longitudinal approach of assessing urban and suburban children's exposure to pyrethroid pesticide. Environ Health Perspect. 2006, 114: 14191423. 10.1289/ehp.9043.View Article
 EPA: EPA's Reanalysis of Key Issues Related to Dioxin Toxicity and Response to NAS Comments (External Review Draft). 2010, Washington, DC: EPA, EPA/600/R10/038A
 Bartell SM, Griffith WC, Faustman EM: Temporal error in biomarker based mean exposure estimates for individuals. J Expo Anal Environ Epidemiol. 2004, 14: 173179. 10.1038/sj.jea.7500311.View Article
 Sherlock JC, Quinn MJ: Underestimation of doseresponse relationship with particular reference to the relationship between the dietary intake of mercury and its concentration in blood. Hum Toxicol. 1988, 7: 129132. 10.1177/096032718800700204.View Article
 Bartell SM, Ponce RA, Sanga RN, Faustman EM: Human variability in mercury toxicokinetics and steady state biomarker ratios. Environ Res. 2000, 84: 127132. 10.1006/enrs.2000.4104.View Article
 Gustafson P: Measurement Error and Misclassification in Statistics and Epidemiology: Impacts and Bayesian Adjustments. 2004, Boca Raton, Florida: Chapman & Hall/CRC
 Ferson S: What Monte Carlo methods cannot do. Hum Ecol Risk Assess. 1996, 2: 9901007. 10.1080/10807039609383659.View Article
 Wedderburn RWM: Quasilikelihood functions, generalized linear models, and the GuassNewton method. Biometrika. 1974, 77: 642648.
 Hardin JW, Hilbe JM: Generalized Estimating Equations. 2003, Boca Raton, Florida: Chapman & Hall/CRC
 Sherlock J, Hislop J, Newton D, Topping G, Whittle K: Elevation of mercury in human blood from controlled chronic ingestion of methylmercury in fish. Hum Toxicol. 1984, 3: 117131. 10.1177/096032718400300205.View Article
 Bartell SM: Statistical Methods for NonSteadyState Exposure Inference Using Biomarkers. PhD Thesis. 2003, University of California, Davis, Graduate Group in Epidemiology
 Liang KY, Zeger SL: Longitudinal data analysis using generalized linear models. Biometrika. 1986, 73: 1322. 10.1093/biomet/73.1.13.View Article
 Prentice RL, Zhao LP: Estimating Equations for Parameters in Means and Covariances of Multivariate Discrete and Continuous Responses. Biometrics. 1991, 47: 825839. 10.2307/2532642.View Article
 McCullagh P, Nelder JA: Generalized Linear Models. 1989, London: Chapman & Hall, 2View Article
 DeGroot M: Probability and Statistics. 1989, Reading, Massachusetts: AddisonWesley, 2
 Stern AH: Estimation of interindividual variability in the onecompartment pharmacokinetic model for methylmercury: implications for the derivation of the reference dose. Regul Toxicol Pharmacol. 1997, 25: 277288. 10.1006/rtph.1997.1105.View Article
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.