# Technical Basis for the Lowest Concentration Minimum Reporting Level (LCMRL) CalculatorOffice of Water (4607M) EPA 815-R-11-001 December 2010 http://water.epa.gov/drink/LCMRL Technical Report December 2010

Page i 

## Table of Contents

1.0 Introduction .......................................................................................................................... 1 1.1 Programmatic Background ............................................................................................... 1 1.2 Scientific Background ...................................................................................................... 2 1.3 Improvements over Previous LCMRL Methodology ...................................................... 3 2.0 General Notation .................................................................................................................. 4 3.0 LCMRL Study Design ......................................................................................................... 5 4.0 Data Conditioning ................................................................................................................ 6 5.0 Estimation of Replicate Variance ........................................................................................ 6 5.1 Introduction ...................................................................................................................... 6 5.2 Modified Hodges-Lehmann Location Estimator.............................................................. 7 5.3 Mean Absolute Deviation Scale Est imator ...................................................................... 7 5.4 Huber Location Estimator and Weights ........................................................................... 7 5.5 Tukey’s Biweight Location Estimator and Weights ........................................................ 8 5.6 Final Location

- Scale Estimator and Weights ................................................................... 9 6.0 Replicate Variance Model .................................................................................................. 10 6.1 D erivation ....................................................................................................................... 10 7.0 Conditional Means Model.................................................................................................. 12 8.0 Conditional Mean Squared Error Model............................................................................ 14 8.1 Estimating Conditional Mean Squared Errors ............................................................... 14 8.2 Fitting Conditional Mean Squared Errors Model ........................................................... 14 9.0 Estimating the LCMRL ...................................................................................................... 15 9.1 Introduction .................................................................................................................... 15 9.2 Prediction Variance ........................................................................................................ 15 9.3 Conditional Distributions ............................................................................................... 16 9.4 Search for LCMRL ......................................................................................................... 17 10.0 Estimating the Modified Hubaux- Vos DL ......................................................................... 20 10.1 Estimating the Critical Level ...................................................................................... 20 10.2 Search for mHV

- DL ................................................................................................... 21 11.0 Summary, Conclusions and Recommendations ................................................................. 22 12.0 References .......................................................................................................................... 23
Appendix A: Computation of MRL ............................................................................................ A-1LCMRL Technical Report December 2010

Page 1 of 24 

## 1.0 Intro duction

1.1 Programmatic Background The Safe Drinking Water Act Amendments of 1996 require EPA to establish criteria for a monitoring program and to publish a list of not more than 30 unregulated contaminants for which public water systems (PWS) are to monitor . The monitoring program will provide a national assessment of the occurrence of these contaminants in public drinking water that will be used to help decide which contaminants may or may not require regulation in the future. In 1999, EPA revised the app roach for unregulated contaminant monitoring in the Unregulated Contaminant Monitoring Regulation (UCMR) (64 FR 50556, USEPA, 1999) and subsequent revisions.
PWSs will be required to monitor for a variety of contaminants under the UCMR. A Minimum Report ing Level (MRL) will be ass igned to each contaminant. Under the UCMR, laboratories will be required to report all occurrences of listed contaminants at concentrations that are equal to or greater than the established MRL.
MRLs represent an estimate of t he lowest concentration of a compound that can be quantitatively measured by members of a group of experienced drinking water laboratories. It is based on a Measurement Quality O bjective (MQO) of 50%-150% recovery of spiked concentrations. Informally, at or above the MRL, competent drinking water laboratories should be expected to obtain 50%-150% recovery or better

. In both the Information Collection Rule (ICR) and the UCMR, the EPA Office of Ground W ater and Drinking Water (OGWDW) specified MRLs and a n accuracy requirement for recovery at the MRL. The MRL was introduced with the new analytes and new methods for implementing the new UCMR.
The Method Detection Limit procedure (MDL; see 40 CFR part 136, Appendix B) was developed by EPA (Glaser et al. , 1981) to establish a decision threshold for detection in lowlevel samples. It is based solely on the standard deviation of repeated measurement of low- level spikes. The MRL differs from the MDL by considering not only the standard deviation of low
conc entration analyses (precision) but also the accuracy of the measurements as they impact achievement of the MQO for spike recovery

. The practical quantitation limit (PQL) was established by EPA’s drinking water program. The PQL is defined as "the lowest concentration of an analyte that can be reliably measured within specified limits of precision and accuracy during routine laboratory operation conditions" (52 FR 25690, USEPA , 1987). The PQL is problematic as a practical matter since at least three diff erent methods have been used to determine PQLs:

1. analysis of performance evaluation (PE) sample data, 2. a risk -based multiple of the MDL, and 3. adoption of Contract Laboratory Program Contract Required Quantitation Limits (CRQLs) based on lowest nonzero standard in a linear calibration.LCMRL Technical Report December 2010

Page 2 of 24 The MRL may be useful as an alternative to the PQL for setting f uture regulatory limits

. A necessary element of determining the MRL is a procedure to estimate the lowest level at which a single laboratory can meet the specified MQOs. EPA has developed a statistical approach for determination of single-laboratory Lowest Concentration MRLs (LCMRLs). This approach us es variance function modeling, regression and modeling instrument response distributions.

The MRL will be determined using a Bayesian bootstrap (BB) (Rubin, 1981) of the LCMRL estimator using the LCMRL study data from each of several experienced drinking water laboratories. The BB replicates that were generated from each laboratory’s data serve to estimate t he distribution of estimated LCMRL values that each laboratory might generate on repeated performance of the LCMRL study. The distribution of pooled BB replicates, generated from the LCMRL study data from a sample of experienced drinking water laboratories, approximates the distribution of estimated LCMRL values which might be generated from the population of experienced drinking water laboratories.
Reasonable values for programmatic required MRLs can be determined based on these simulated distributions. Drinking water laboratories will confirm that they are capable of meeting a required MRL during their Initial Demonstration of Capability (IDC) for an analytical method

. A technical description of the MRL procedure is given in Appendix A.



### 1.2 Scientific Ba ckground The Lowest Concentration Minimum Reporting Level (LCMRL) is defined as the lowest spiking concentration such that the probability of spike recovery in the 50% to 150% range is at least 99%. The algorithms and procedure presented in this report provide a means of reliably estimating the LCMRL.
The LCMRL calculator constructs conditional (on spiking concentration) mean and variance models for analytical measurements. For analytical methods which cannot give negative ‘raw’ measurements or negative results (except possibly at low level due to a negative intercept in the calibration function), the gamma distribution is used as the distribution of measurements

. Reported negative measurements due to a calibration artifact are set to zero. Typically, t hese analytical methods are chromatographic methods in which the response is integrated valley to valley.

Other classes of methods can obtain ‘raw’ measurements at very low concentration that are negative. These methods may be spectroscopic or radiochemical methods and often involve subtraction of ‘dark current,’ reference cell readings or background

. In these cases the normal distribution is used as the distribution of measurements.
The conditional mean and variance models then specify the parameters for the conditional measurement distribution as a function of spiking level. It is this model of the distribution of repeated measurements at a given spiking level that is used to estimate the LCMRL.LCMRL Technical Report December 2010

Page 3 of 24 The critical level , L c is defined as an upper percentile of the distribution of repeated measurements at zero true concentration (see Section 10.1, and Currie, 1968 ). The modified Hubaux- Vos Detection Limit (mHV -DL) is defined as the spiking concentration that has a 95 % probability of exceeding the 95% L c . The conditional distributional model constructed for determining the LCMRL is also used to determine the mHV- DL.



### 1.3 Improvements over Previous LCMRL Methodology The currently proposed LCMRL methodology offers several improvements over the previous LCMRL (USEPA, 2004).
The first improvement is the fundamental use of a probability-based definition for the LCMRL in deriving the estimation algorithms . As described above, the LCMRL is defined as the lowest spiking concentration such that the probability of spike recovery in the 50% to 150% range is at least 99%. The previous LCMRL calculator used a graphical technique which did no t completely agree with the probability - or MQO -based definition of the LCMRL. The probability-based definition leads to an approach based on characterizing the distribution of the analytical response as a function of the spiking level.
The second improvement is that the piecewise linear replicate variance function in the first LCMRL calculator has been replace d by a smooth param etric function. The piecewise linear model was based on linear interpolation of pointwise sample variances. The proposed methodology replaces the ordinary sample variances with robust estimates of variance based on modern statistical methods (M

- estimatio n). This enables the potential impact of outliers to be effectively eliminated without requiring outlier testing and removal.

The thir d improvement is that t he robust replicate variance estimates are fit to a model based on an additive and multiplicativ e error structure for trace level analytical chemistry measurements.

The fourth improvement is that the average response model (regression of measured onto spiked concentrations) is fit using iteratively re weighted least squares ( IR LS) in which the weights are determined not only by the reciprocal of the robust variances but also by M- estimation weights determined as part of the robust variance estimation process. This again enables the potential impact of outliers to be effectively eliminated without r equiring outlier testing and removal.
The fifth improvement is in considering lack of fit for the average response model as a component of prediction variance. This is done by computing robust estimates of mean squared error (MSE) and fitting a parametr ic model for conditional MSE, similar to what is done for replicate variance. The prediction variance function is constructed using the pointwise maximum of the replicate variance and conditional MSE functions. This is multiplied by a function that inclu des the additional uncertainty inherent in the estimation of the average response model. This is a conservative approach based on the author’s perception that the greater danger arises from underestimating the prediction variance.
The sixth improvement is that the gamma distribution is used as a response model for analytica l methods which do not give negative results at low level. The normal distribution is used as theLCMRL Technical Report December 2010

Page 4 of 24 response model for analytical methods which can give negative results at low level. The degrees of freedom in the normal case are estimated from the approximate degrees of freedom for the variance predictions given by the replicate variance model.
The seventh improvement is using the conditional response distribution estimated to determ ine the LCMRL to estimate L c and the detection limit using a modification of the Hubaux- Vos approach.



## 2.0 General Notation Below are definitions of key variable notations used throughout this paper. Additional notations included in this document are described as they are introduced.

m is the number of spiking levels.

n is the total number of observations.

n i is the number of replicate observations at the i th spiking level.

x i is the spiking concentration for the i th spiking level.
y ij is the measured concentration for the j th measurement at the i th spiking level.
Y

ij is the j th measurement at the i th spiking level represented as a random variable.

w ij is the weight for the j th measurement at the i th spiking level.
y is the n x 1 vector of measured concentrations.

p is the number of model parameters. **β** is a p x 1 vector of regression coefficients (model parameters).

E (.) is the expected value function for a random variable. var(.) is the variance function for a rand om variable.

i **μ** is the expected value of Y ij .

2 i **σ** is the variance of Y ij . Z is a design matrix.

i Z is a design vector (row of the design matrix) for the measurements at the i th spiking level

. CV is coefficient of variation (also known as relative standard deviation).
L c is the critical level

.

C **γ** is the coverage probability for the critical level , L c .

D **γ** is the coverage probability for the m HV

- DL.

Q **γ** is the coverage probability for the quality specification interval (50%- 150% re covery) at the LCMRL

.LCMRL Technical Report December 2010

Page 5 of 24 

## 3.0 LCMRL Study Design Because the LCMRL methodology is a regression-based methodology, LCMRL study designs involve replicate spiking at multiple levels. The relevant design parameters are:

**•** the number of spiking levels,

**•** the spacing of spiking levels,

**•** the number of replicates at each spiking level (not necessarily equal) , and

**•** inclusion of method blanks.

The relevant issues ass ociated with the design of detection and quantitation studies include the impact of the study design on:

**•** identification of nonlinearity and the correct model form for the average response function,

**•** parameter estimation in the replicate variance and condit ional MSE model, and

**•** parameter estimation in the average response function model.
Because the LCMRL methodology is complex and multistage, definitive recommendations for study designs are not available at present. This is an issue for further development.
Hubaux and Vos (1970) study a variety of standards spacings, which they refer to as repartitions of standards. They study linear spacing, parabolic spacing, two point spacings (that is, two spiking levels) and three point spacings. They recommend th at the ratio of the maximum to minimum (nonzero) spiking level be about 10. They also recommend having all the replicates (save one) at the lowest and highest spiking levels. This is the ir three point design. This report rejects that particular recommen dation for the following reasons.
Both the average response and replicate variance models are regression models. In design for regression, having more levels of the predictor (the spiking level in this case) is much more important than having more replication at each level. However, for the LCMRL, estimation of the replicate variance function is also vital. Therefore, some replication is also essential. Under Gaussian sampling, the coefficient of variation of the ordinary sample variance is 1 for three replicates , 82.0 3 2 **≈** for four replicates , 63.0 5 2 **≈** for six replicates, 53.0 7 2 **≈** for eight replic ates and 47.0 9 2 **≈** for ten r eplicates.
Nonlinearity of the response function at low concentration is a well known phenomenon. I t is very important to characterize nonlinearity at low level in order to accurately estimate the LCMRL and DL because the mean response function at low level is a key component of the LCMRL/mHV -DL methodology. For this reason , we allow a polynomial model, up to a cubic, for the mean response. In order for the coefficients of a cubic polynomial to be identifiable (estimable), there must be at least four spiking levels.
It is also very important to estimate the variance functi on as accurately at low level as possible. For these reason s, this report recommends the parabolic design for spiking levels with as many levels as feasible and at least four replications per spiking level.LCMRL Technical Report December 2010

Page 6 of 24 The spiking levels may be computed as follows 

## ( ) 

# ( ) 2 1 1 1 1 

### , for 1, , , i i m m 

### xx x x i m

**−**

**−** 

### =+ − **=**  (1)

where x 1 and x m are the designed minimum and maximum spiking levels.
The proposed LCMRL calculator recommends seven spiking levels with four replicates per level for a total of 28 analyses. Fitting a cubic polynomial in the mean response model requires a minimum of four spiking levels. The minimum values that the LCMRL calculator will accept are four spiking levels and three replicates per level. Given that the size of the LCMRL/DL study must be limited by cost and resource availability, a tradeoff must be made between estimating replicate variance and MSE at individual spiking levels and estimating the regre ssion models for mean response, replicate variance and MSE. The recommended numbers of spiking levels and replicates are a reasonable compromise.



## 4.0 Data Conditioning Because the LCMRL methodology is a complex, multistage , model-based design, care must be taken in screening the input data for use in the model development. The model estimation procedure issensitive to departures from a monotonic variance model. Having non-zero spiking levels with nonzero replicate variances is problematic for modeling. Therefore, data values that would contribute to invalidating these assumptions are dropped from the estimation procedure. Data are dropped from the estimation procedure for non-zero spiking levels where 50% or greater of replicate responses are zero. If less than 50% of the replicate responses of a non-zero spiking levels are zero, then the zero responses are replaced by the minimum non-zero replicate response within the spiking level or within a lower spiking level.
After data conditioning the procedure determines if the requirement of a minimum of four spiking levels and three replicates is achieved before estimation of the LCMRL.



## 5.0 Estimation of Replicate Variance

5.1 Introduction Because of the potential for outliers, r eplicate variance at each spiking level is estimated using weighted variances using weights computed by M- estimation . As a byproduct of this process, robust weights are available for estimating the conditional mean (regression) model. Since at each spiking level the number of results may vary from 3 to perhaps 10 or more, this is a very challenging undertaking.
The objective is not to ‘identify’ outliers as such but merely to down- weight them sufficiently that they do not overly influence the means model, replicate variance model or conditio nal MSE model, without overly compromising efficiency. In contrast, testing and removal of outliers constitutes a data weighting scheme where the only possible weights are 0 or 1 n , where n is the number of retained observations. This is very inflexible compared to weighting based on M-LCMRL Technical Report December 2010

Page 7 of 24 estimation. Not only is testing and removal of outliers awkward procedurally, it causes a loss of statistical efficiency.
A modified Hodges- Lehmann estimator is used for the initial resistant location estimate

. The scaled mean absolute deviation from the Hodges

- Lehmann estimator is used as the resistant scale estimator. Using this initial location

- scale estimation pair, the Huber M- estimate of location and its associated weights are computed. A robust scale estimator is computed based on the weighted mean absolute deviation from the Huber location estimator and using the Huber weights. The Huber location- scale estimates are then used as an initial location and robust scale estimates for Tukey’s biweight procedure.
Finally, the Tukey’s biweight weights are used to compute location and variance estimators at each spiking level. The weights are also later used in computing the conditional mean response (regression) model.



### 5.2 Modified Hodges

- Lehmann Location Estimator The Hodges

- Lehmann location estimator (Randles and Wolfe, 1991) is the median of the set of means of all pairs of observations. The modified Hodges-Lehmann location estimator (modified for this application) at the i th spiking concentra tion combines the median of the measurement data with the set of means of pairs of data and takes the median of the resulting set of numbers. For a set of n i observations 

# { } 1 , , i i in y

 y for the i th spiking level, this is given by 

# { } , 1 median ,median , , 2 i ij ik HL i i in j k y y m y y **< +**  ** **

**=  **  ** **  

# 

. (2)



### 5.3 Mean Absolute Deviation Scale Estimator This estimator is based on the long known (Kelley, 1921) Absolute Deviation (AD) scale estimator . It uses the mean absolute deviation from the Hodges- Lehmann estimator (instead of from the sample mean ). It is computed as

1

, 1 1.4826 i n i ij HL i n j AD y m

**= =⋅ −**  

## ∑

. (3)

The factor 1.486 is a normalization constant, used so that the expected value of the AD in normal samples equals the population standard deviation.



### 5.4 Huber Location Estimator and Weights The Huber estimator (Hoaglin e t al. , 1983) is calculated iteratively as a weighted estimator ( w
estimator) using a relative convergence criterion of 10 -6 . The process is started with a resistant initial location estimate (0)

, H i T given by the modified Hodges-Lehmann estimator and a resistant auxiliary scale estimate given by the modified AD estimator, i AD . This is a departure from theLCMRL Technical Report December 2010

Page 8 of 24 usual practice, in which the median is used as the initial location estimate and the median absolute deviation ( MAD ) scale estimator is used as the auxiliary scale estimator for the Huber location estimator. The modification is prompted by the very small sizes of the data sets in this application.
The Huber weights and estimator are determined by a tuning constant c H according to the equations below. The value of the parameter c H is taken to be 1, which is a standard value often used in practice. The index k tracks the iteration number. For a set of n i observations 

# { } 1 , , i i in y y

 at the i th spiking level the iteration proceeds as follows

**( ) ( )** 

**( ) ( )**  **( )**

**( )** 

# ( **)** 1

,

,

, ,

, 1 1 , 1, ,
1, if , 1, ,

, if k ij H i

i

i i y T k ij i AD

k ij H H k ij i k k H ij ij H

n n k H k H k H i ij ij ij j j u jn

u c w j n c u u c T wy w **−**

**−**

**= =**  **= = **  **≤ ** 

**= =** 

** >**  ** ** 

**=** 

## ∑ ∑ 

 (4)

until the estimate **( )**

, k H i 

### T converges, at whi ch time the weights are assumed to have converged based on a relative convergence criteria of 10 -4

. The Huber M- estimate of location is denoted by

, H i T and the associated weights by H ij w

.

The Huber scale estima tor, 2

, H i s (equation 5), is computed as the weighted variance around the Huber location estimate using the Huber weights. This formula is derived from the variance of a weighted mean. It also follows from an approximation to the A- es timators of scale (Lax, 1985, Section 4.6.2) based on the asymptotic variance of M- estimators (Hoaglin et al. , 1983). This completes the circle, since M

- estimators can always be represented as weighted means ( Westimators). 

# ( ) **( )**

2

, 1

2 2

, , 1

, 1 i

i n H w i i ij j

n H i H i ij ij H i j w i nn w

n s wyT n **=**

**=  ** 

**= −  **  ** **

**= −**  

## ∑ **∑** (5)



### 5.5 Tukey’s Biweight Location Estimator and Weights The Tukey’s biweight M- estimator (Hoaglin et al., 1983; Horn , 1988 ) is calculated iterativ ely as a weighted estimator using a relative convergence criterion of 10 -4

. The biweight is more robust statistically in certain respects than the Huber estimator. It is very efficient and can down -weight extreme observations to zero. On the other hand, it is not as robust computationally as the Huber estimator. Without a good auxiliary scale estimate and aLCMRL Technical Report December 2010

Page 9 of 24 good starting location estimate, it can misbehave. To address this issue t he process is started with a robust initial location estimate T (0) BW given by the Huber estimator and a robus t auxiliary scale estimate given by Huber weighted AD scale estimator, s H i , . This is a departure from t he usual practice, in which MAD is used as the auxiliary scale estimator . It is prompted by the very small sizes of the data sets in this application.

The Tukey weights are determined by a tuning constant c BW according to the equations below. The value of the tuning parameter c BW is taken to be 9. The index k tracks the i teration number. For a set of n i observations 

# { y i 1 ,  , y in i 

# } at the i th spiking level the iteration proceeds as follows

**( ) ( )** 

**( )** 

# ( ) 

# ( ) **( )**

**( )**

**( )** 

# ( **)** 1

,

, 2 2

,

, ,

, 1 1 , 1, ,

1 , if 1 , 1, , 0, if 1 k ij BW i
BW H i

i i y T k ij i c s

k k ij ij BW k ij i k ij

n n k BW k BW k BW i ij ij ij i i u jn

u u w j n

u T wy w **−**

**−**

**= =**  **= = **  **− < ** 

**= =** 

** **  **≥ ** 

**=** 

## ∑ ∑ 

 (6)

until the estimate 

### T **(** k **)** BW , i converges, at which time the weights are assumed to have converged. The location estimate is denoted by T BW , i and the weights by w BW ij . Note t hat the definitions of

u **(** k **)**

ij are similar but different in the two cases (Huber and biweight).



### 5.6 Final Location

- Scale Estimator and Weights The final robust weights are the Tukey biweights. The robust location and variance estimates for the i th spiking level are computed as a weighted mean , m w i , , and weighted variance, s 2 w i

, , of the data as indicated below where nw i is the degrees of freedom associated with the weighted observations. 

# ( ) 

# ( ) **( )** 2

, 1

, 1 1

2 2

, , 1

, 1 i

i

i n BW w i i ij j

n n BW BW w i ij ij ij j j

n BW i w i ij ij w i j w i nn w

m wy w

n s w ym n **=**

**= =** 

**=  ** 

**= −  **  ** **

**=**

**= −**  

# ∑ **∑ ∑ ∑**  (7)

As noted above, th is is the formula for a weighted variance , derived from the variance of a weighted mean. It is also an approximation to the A

- estimators of scale based on the asymptotic variance of Tukey’s biweight location estimator

.LCMRL Technical Report December 2010

Page 10 of 24 C haracteristics of the Huber and biweight estimators are combined by using the Huber estimates as scale and st arting location in computing the biweight. The reason for doing this is that for M

estimators with monotone ψ -functions, like the Huber, the equation defining the M- estimator is guaranteed to have a unique solution. For redescending estimators, like Tukey’s biweight, uniqueness is not guaranteed. This problem is particularly worr isome given the small data sets being used (results at each spiking level). Since the computer code for LCMRL/mHV- DL estimation must run unsupervised, the Huber estimate is a safe starting point for the biweight iterations.



## 6.0 Replicate Variance Model

6.1 Deriv ation The replicate variance model, which models replicate variance as a function of spiking concentration, is based in part on the paper by Rocke and Lorenzato (1995). This paper identifies sources of both additive and multiplicative random error in trac e level analytical chemistry measurements. Corresponding to this conceptual model, the LCMRL procedure fits constant plus power models to replicate variance and conditional MSE model s. The constant term estimates the additive sources of variance and the power model term estimates the multiplicative sources. The power model term accommodates error variance behavior in higher concentration ranges that varies from near constant variance to shot noise to constant coefficient of variation (CV).
The variance at extremely low concentrations may be distorted by processes such as thresholding, smoothing, rounding and truncation. Smoothing, which is often used in spectroscopic, chromatographic, and radiologic methods, reduces measurement variance, especially fo r low

- level measurements. Rounding and truncation also reduce measurement variance.
Response thresholding, which is used in chromatographic methods to manage data volume, can either reduce or increase measurement variance depending on the situation. Thresholding effectively replaces noise (except for the extreme upper tail of the noise distribution) and some low level signal with zeros. Thresholding then functions as a Type II censoring mechanism with an unknown censoring point (unless the laboratory reports the instrument settings)

. Instrument thresholds have both direct and indirect impact s on estimating detection and quantitation limits. The main direct impact is that it is not possible to estimate the standard deviation of measurements at zero. However, by definition, the standard deviation at zero is required to calculate L c (Currie, 1968). The EPA MDL procedure (Glaser et al., 1981) was constructed to deal with this problem by providing a way to estimate a standard deviation at a low concentration, and includ ed instructions for determination of a concentration as close to zero as is possible that will generate a measurement.
Theoretically, in the absence of such distortion, low- level measurement variance (even for method blanks) would always be positive. For this reason, method blank data are not used in estimating the replicate variance model.LCMRL Technical Report December 2010

Page 11 of 24 The replicate variance model takes the form 

## ( ) 2

, v c
v v x a bx **σ = +**  (8)

where x is the spiking concentration.
The replicate variance model parameters are estimated using a constrained iterated downhill Simplex Method minimization routine (Nelder and Mead, 1965) with starting values for b and c found by first fitting the robust variances using log- scale simple regression. Because of the potential for outliers, the ordinary sample variance is not used as the response in the fitting

. Instead a robust estimate of variance is compute d bas ed on M- estimation , as described above.
The starting values 0 b and 0 c , for b and c, are first found by fitting a reduced variance model of the form 

## ( ) 2 ln[ ] ln( ) ln( ),
v v x bcx **σ = +** 
using ordinary regression using the robust variances 2

, w i s at each spiking level i x

. Following this, the full replicate variance model is fit using a n iterated constrained downhill Simplex Method with starting values defined as



## ( ) 

## ( **)** 2 2

, , 1 2

,

1

0

0 max ,1 8

min max 0, ,2 wi wi i
s w i i

s

s n s a e n b b

c c **=**

**=  **  ** **
**= −**  ** **

**=**

**=** 

## ∑ **∑**

a nd a constrained loss function defined a t each iteration j of the optimization as



# ( )
**( )** 1

12

2 2 2

, 2

2 th where 10 , if 0 or 0 or 0 or 2

, otherwise
, at the iteration of the optimization. j n i i

jjjj

i i i w i i

c k kk loss e

abcc

e s n

x a bx k **σ**

**σ**

**σ =**  **=**
** <<<>**  ** ** 

**= ** 

**− **  ** ** 

**= +**



## ∑

T he optimization continues for estimates of a, b and c until the loss function is less than 10 - 16

. Normal variances (that is, the statistic computed from data) have their sample variances proportional to their means. The robust variances used in the model behave similarly. It is thisLCMRL Technical Report December 2010

Page 12 of 24 that motivates the choice of the functional form of the loss function. The large loss incurred when parameter estimates are out of valid ranges forces the Simplex Method to constrain the parameter values to valid ranges.



## 7.0 Conditional Mean s Model The conditional means model is estimated using a robust estimation procedure that implements a nested Iterative Reweighted Least Squares (IRLS) procedure. As part and as a consequence of the iterations , the conditional mean squared error model and the robust weights are also updated until convergence of the procedure . The IRLS procedure can handle non- st andard regression assumptions such as non-homogenous variances and non- normal errors . This method is also resista nt to outliers in the data set by using weights that are derived using Tukey’s bis quare method (See Section 4.5).
The general conditional mean s model is 

## ( ) 

# ( ) 

# ( ) 2 3 E

1, , , T i ij i

T i ii i xY Z
Z xx x **μ β = =** 

**=** (9)

To solve this in the IRLS framework, it is necessary to solve

1 1 ˆ ( )( ) ( ) , n n i i i i i ii i i w u y z z w u ez **β**

**= =**  **− =** 

# ∑ ∑

w here e i is the residual at each point and w ( u i ) are the Tukey bis quare weights (estimated independently here for the means model). The Tukey weights, defined in Section 5.5, are in this
case a function of the x i , the predicted values as an estimate location , the conditional mean squared error as an estimate of scale and a tuning constant of 9. The conditional mean square error is used instead of the traditional variance to allow for the model lack of fit.

I t is fir st necessary to have in itial starting values f or 

### β ˆ and the conditional mean square error (s ee Section 8.0). Using th e in i tial ˆ **β**
, say 0 ˆ **β**
, the i nitial predicated values, residuals and the conditional mean square error (s ee Section 8.0) are calculated

. The initial 0 

### ˆ **β** is found by performing weighted least squares regression using the robust weight matrix W (see Section 5.6). The solution for 0 ˆ **β** using weighted least squares regression i s given below.LCMRL Technical Report December 2010

Page 13 of 24 

# ( ) 

# { **} (**  **)**

1 1 0

1,1 ,

1,1 1,2 1, ,1 ,

1

1

1 ˆ diag , ,

, ,, ,, ,,

rows

rows m

m T

m n
T

n m mn

T

T

T

m

m T

m y

w w

y yy y y y
Z

n Z

Z

n Z **β −**  **=**

**=**

**=  **  ** **  **  **  **  **  **  **  ** **

**=  ** 

**  **  **  **  **  **  **  **  **  ** 

ZWZ ZW
W

Z 









The nested IRLS is implemented for estimating 

### ˆ **β** as follows. Parameters epsilon and maxIter are 10 -6 and 100, respectively.

1) From the initia l 2) For the outer loop from the initial residuals compute the initial estimate of the conditional MSE function, 3) For the inner loop compute the initial Tukey weights W 0 using the predicted values and cond itional MSE function 4) Us e weighted least squares to obtain the estimate 5) Us e the parameter estimate 6) G o b ack to step 3 until the estimate 0 ˆ **β** calculate the predicated values and residuals.

0 mse

.

0 mse and a tuning constant of 9.

ˆ i **β**

.



### ˆ i 

### β to obtain new predicated values and residuals.

1 ˆ ˆ | | i i **β β**

**− −**  converges pointwise such that the

max imum pointwise difference is les s than epsilon or the number of iterations exceeds maxIter . S et 

### ˆ ˆ i 

### β β **=** and go to step 7. 7) Go to step 2 and compute a new conditional mean square error function mse j using ˆ **β** and the residuals from step 5. 8) I f the maximum pointwise difference of 0 

### ˆ ˆ | | **β β**

**−** (over the outer loop) is greater than the convergence tolerance and the maximum number of outer iterations has not been exceeded, set 0 ˆ ˆ **β β =**  and 0 j mse mse **=** and go back to step 1. Otherwise , stop and return
the current values 

### ˆ **β** and j mse

.
Final m odel selection (linear, quadratic or cubic) is accomplished using Mallow’s C p (Draper and Smith, 1981). Since the constant term in the model may be negative, it may be possible for mean concentrations near zero to have a negative mean estimate. Th erefore , the mean response function at a spiking concentration, x , is defined asLCMRL Technical Report December 2010

Page 14 of 24 ˆ 0, if 0 ( )

ˆ ˆ , if >0 z x z z **β μ**  **β β **  **≤ ** 

**= **  ** ** 



## 8.0 Conditional Mean Squared Error Model

8.1 Estimating Conditional Mean Squared Errors Since the sampling distribution for repeated measurements at a spiking level is based on the conditional means (regression) model and variance about the conditional mean, the contribution to variance by lack of fit must be accounted for. To account for the effect of lack of fit for the regression model, a model is fit for conditional mean squared error (cMSE) as a function of spiking level.

First , robust estimates are made of MSE at each spiking level. T his process is anal ogous to the estimation of replicate variance except that the raw residuals from the conditiona l means model at each spiking level are used as the data . The robust c MSE for the i th spiking level, x i , is calculated as

2 2

, , i wi wi cMSE m s **= +** (10)
with the quantities calculated as described in S ection 5.6 but using the regression residuals. Clearly, the c MSE should be larger than the replicate variance at each spiking level.



### 8.2 Fitting Conditional Mean Squared Errors Model The conditional mean squared error model, which models mean squared error as a function of spiking concentration, is assumed to be in the class of constant + power models

2

, e c i e ei a bx **τ = +**  (11)

where i = 1, 2, ..., m indexes the spiking levels, x i . The cMSE model is fit in exactly the same manner as the replicate variance model. Because of the model fitting process however, there is no guara ntee that the expected value from the cMSE model is larger than the expected value of the replicate variance model at every spiking level.LCMRL Technical Report December 2010

Page 15 of 24 

## 9.0 Estimating the LCMRL

9.1 Introduction To construct an estimated sampling distribution for repeated measurements at an arbitrary spiking concentration, a mean function, a variance function and a parametric distributional family are required

. The mean function is estimated as described in S ection 7.0.
As described in Section 1.2, the gamma distribution is used for analyt ical methods that give only
non- negative measurements , and the normal distribution is used for methods that can give both positive and negative measurements at low concentration. The gamma distribution is appropriate in the first case since it is the maxi mum entropy distribution (Kagan et al. , 1973) in the class of continuous distributions which support only non- negative values and have a specified mean and variance (or specified mean and geometric mean). The normal distribution is the continuous maximum entropy distribution with specified mean and variance and support on the real line.



### 9.2 Prediction Variance As described in S ection s 6.0 and 8.0, models are constructed for both replicate variance and conditional MSE (of residuals from the means model) as functions of spiking concentration. For estimation of the LCMRL and HV-DL, what is really needed is a model for prediction variance as a function of spiking concentration. Three additional considerations allow construction of a robust yet conservative model for prediction variance.

The first is that , as discussed in S ection 6.1, the variance at extremely low concentrations may be distorted by processes such as thresholding, smoothing, rounding and truncation. Therefore, a conservative prediction variance function should never give a zero prediction variance, not even at zero spiking concentration.
Because the constant in either the replicate variance (equation 8) or cMSE (equation 11) models may be zero, a parameter called minVar has been created . The minVar parameter is set to the constant a if it is nonzero. Otherwise it is set to the mean of the robust variance s at the two lowest nonzero spiking levels. The minVar parameter is denoted by 2 min **σ** and 2 min **τ** fo r the replicate variance and cMSE models, respectively

. The second consideration is that as noted in S ection 8.1 , although the replicate variance at each spiking level should always be smaller than the c MSE, it is possible for the replicate variance model estimate to be larger at some spiking concentrations than the cMSE model estimate. A conservative approach is to use the point- wise maximum of these functions.
The third consideration is that the prediction variance is associated in a sense to the reg ression (conditional means model). The prediction variance should be augmented with regression model uncertainty just as it is in ordinary least squares (OLS) regression (Draper and Smith, 1981). This makes intervals calculated using the prediction varia nce comparable to tolerance intervals.
These considerations guide the construction of the prediction variance function:LCMRL Technical Report December 2010

Page 16 of 24 

## ( ) **( )** 

# ( ) **( ) ( )**  

# ( )
**( ) ( ) ( )**  

# ( **) ( )**  **( ) ( ) ( )**  

# ( **) ( )**  1 1
2

2
1
2

2

2 2 2 min min

2 2 2 22 2 1 1

max ,

max , , if max , 0

max , , , otherwise v

e m j j n j

m j j j

c
v v

c
e e

mm

v e pred

mm x nx

x x f x n nx x x a bx
x a bx

x x fx a a x x x fx **σ**

**τ**

**σ στ**

**σ τ σ**  **σσ τ =** 

**= =** 

**− =+ +** 

**−**

**= +**

**= +**

**= **  **⋅ > ** 

**= ** 

**⋅ **  **** 

## ∑
**∑**

(12)

Prediction variance estimate s from this function are assumed to have p d nmn

**=− −** degrees of freedom, where n p is the maximum of the number of parameters fit in the variance model and the cMSE model.



### 9.3 Conditional Distributions Maximum entropy is one method of characterizing probability distributions (Kagan et al. , 1973). Often there is a unique distribution that maximizes entropy while satisfying conditions of allow able values for a random variable and constraints on specified moments (or on other functionals). This unique distribution satisfies the given conditions while implying absolutely no additional information. Maximum entropy distributions are very important when working with scientific problems that have physical cons traints. They can honor those constraints without adding any additional unwarranted information.
The gamma distribution is the maximum entropy distribution in the class of continuous distributions with support (that is, allowed values) on the nonnegativ e real line and specified mean and geometric mean (Kagan et al. , 1973). This meshes ideally with the model for lowlevel analytical error variance described by Rocke and Lorenzato (1995) and discussed in S ection 6.1. It has both additive and multiplicative error components, corresponding to the linear and log- scale moment constraints in the gamma distribution.
It can easily be shown (Kagan et al. , 1973) that in this case, the constraint on mean and geometric mean is equivalent to a constraint on mean a nd variance. Thus, for specified conditional mean and variance functions with posi tive range the principle of Maximum Entropy generate s a family of gamma distributions. Therefore, the gamma distribution is used as the distributional model for response for methods which do not give negative results.
The normal distribution is the continuous maximum entropy distribution with specified mean and variance and support on the real line (that is, negative values are allowed). The normal distribution is therefore used as the distributional model for response for methods which can give negative results. In the case of the normal model for repeated measurements, the conditionalLCMRL Technical Report December 2010

Page 17 of 24 mean and prediction variance functions directly give the parameters for the repeated
me asurement model. Since the variance function estimates have an associated degrees of freedom d that will typically be less than 30, p robability calculations are made using the t

distribution with d degrees of freedom.
In the case of the gamma distributi on (used for methods which cannot give negative results), the conditional gamma distribution parameters are computed as 

### ( ) ( ) ( ) **( ) ( ) ( )** 2 2

2

. pred

pred xx x

x xx **α μσ βσ μ**  **=**

**=** (13)



### 9.4 Search for LCMRL An iterative search is conducted for the spiking level which satisfies the definition of the LCMRL. In the case of the normal distribution model, the LCMRL is the x which satisfies



### ( ) ( ) 

# ( **) ( )**  **( ) ( )**  **( )** 2 1.5 0.5 Pr 0.5 1.5 , pred d d Q pred pred xx xx xY x x x T T x x **μ μ μ σ**  **γ σ σ**  ** −  −  ≤ <**  **= −**  **=   **  **  ** (14)

p d nmn

**= − −** 

wher e **( )** d T

**⋅** is the cumulative distribution function for the t

- distribution with d degrees of freedom and 0.99 Q **γ =**  .

I n the case of the gamma distribution model, the LCMRL is the x which satisfies


### ( ) ( ) 

## ( **) ( ) ( )**  

## ( **) ( ) ( )**  

## ( **)** Pr 0.5 1.5 , 1.5 , 0.5 , Q

xY x x x
F xx x F xx x **α β αβ αβ γ**  **≤ < =** 

**− =** (15)

where 

## ( )

, F **α β**

**⋅** is the gamma cumulative distribution function with parameters α and β

.

F igure 9 -1 below shows an example of the coverage probability for the MQO recovery interval (50%-150%) as a function of spiking level. Figure 9 -2 below i s an example of an LCMRL plot showing the raw data, the estimated average response curve, th e recovery limits, a 99% prediction envelope for the response and the LCMRL. The text on the figures also gives the LCMRL. These figures were created using the LCMRL calculator software application produced for laboratory use.

I n this example, the spiki ng range is relatively narrow and the replicate variance and conditional mean squared error are relatively constant across the range. Figure 9 -3 below shows an example of an LCMRL plot for an analyte with error (replicate variance and conditional MSE) inc reasing with spiking concentration.LCMRL Technical Report December 2010

Page 18 of 24 Figure 9-1 : Example Coverage Probability for MQO Recovery Limits

Figure 9-2 : Example LCMRL Plot w ith Constant Error Varianc e

0 0.05 0.1 0.15 0.2 0.25 0 0.05

0.1 0.15

0.2 0.25

0.3 0.35

0.4 1,4-dioxane--LCMRL Plot

Spike Concentration ug/LMeasured Concentration ug/L Data LCMRL = 0.042 ug/L Y=X Regression 50-150% Recovery Lower/Upper Prediction Limits

Spik ed Concentration ug/L 0 0.1 0.2 0.92 0.94 0.96 0.98 1 LCMRL = 0.042 ug/L Qual Lim: 50-150% Coverage Prob: 0.99 1,4-dioxane--QC Interval Coverage Plot

Spike Concentration ug/LProbability of QC Interval Coverage

Spiked Concentration ug/LProbability of QC Interval CoverageLCMRL Technical Report December 2010

Page 19 of 24 Figure 9-3 : Example LCMRL Plot with Increasing Error Variance

Figure 9 -4 below shows the increasing MSE as a function of spiking level. The points a re the robust estimates of MSE. The line is an approximation to the constant + power model fitted to the cMSE. 0 100 200 300 400 500 600 700 0 200 400 600 800 1000 Trichloroethene--LCMRL Plot

Spike Concentration ug/LMeasured Concentration ug/L Data LCMRL = 57 ug/L Y=X Regression 50-150% Recovery Lower/Upper Prediction Limits

Spiked Concentration n g/L Trichloroethylene

-- LCMRL Plot

n g/LMeasured Concentration ng/L

mHV

- DL = 18 ng/L L c = 19 ng/L Equation: y = 8.9221 + 1.0564 * xLCMRL Technical Report December 2010

Page 20 of 24 Spiked Concentration ng/L Figure 9-4 : Conditional MSE and fitted model showing increasing error with concentration



## 10.0 E stimating the M odified Hubaux

- Vos DL The proposed methodology differs from and improves on the original Hubaux-Vos methodology (Hubaux and Vos, 1970). First, the proposed method accommodates nonconstant variance in a very robust manner. Hubaux and Vos were aware of this consideration but did not include it in their methodology. Secondly, in cases w here the instrument response cannot be negative, the proposed methodology uses the gamma distribution (which is the Maximum Entropy distribution in this situat ion) to model the response distribution.



### 10.1 Estimating the Critical Level L c was defined by Currie (1968) as an upper percentile of the distribution of repeated measurements at zero concentration. For Currie , L c represents an a posteriori decision limit f or detection. As noted in S ection 6.1, although in many situations it may be impossible to sample this distribution directly in an accurate way , it is possible to estimate using regression methods. This was the general approach of Hubaux and Vos (1970).LCMRL Technical Report December 2010

Page 21 of 24 As is the case for the LCMRL ( S ection 9.1), the methodology used to estimate L c in the mHV

DL method distinguishes two cases. In the first case, in which blank and low- level measurements can give negative instrument responses, the normal distribution is used to model the response distribution. Then L c is estimated as 

## ( ) **( )** 

## ( )

0 0 C pred d C

p L t d nmn **μσ γ =+ ⋅** 

**=− −** (16)

where 

## ( ) d C t **γ** is t he C **γ** -quantile of the t

- distribution with d degrees of freedom (see S ection 9.2), and 0.95 C **γ =** 
.
In the second case , in which blank and low- level measurements cannot give negative instrument responses, the gamma distribution is used to model the response distribution. Then , assuming that 

## ( ) 0 0 **μ >**  , L c is estimated as 

## ( ) **( )** 

# ( **)** 1 0, 0 C C L F **γα β**

**− =**  (17)

where 

## ( ) 1

, F **α β**

**−**
**⋅** is the quantile function of the gamma distribution with parameters α and β

.
However, in the case that 

## ( ) 0 0 **μ =**  , the gamma distribution is not tenable. In this case, for want of a better alternative, a half

- t

distribution is used. Then L c is estimated as 

# ( ) **( )** 1 2 1 0 C C d pred L t **γ σ**  **− =−⋅**  (18)

where 

## ( ) d t **⋅** is the quantile function of the t distribution with d degrees of freedom

.



### 10.2 Search for mHV

- DL An iterative search is cond ucted for the spiking level which satisfies the definition of the mHV

DL. In the case of the normal distribution model, the mHV- DL is the x which satisfies



### ( ) ( ) 

# ( **) ( )**  **( )** 2 Pr , C C pred d D pred L x YL x x T x **μ μ σ**  **γ σ**  **−  ≥ ==**  **   **  (19)

where **( )** d T

**⋅** is the cumulative distribution function for the t

- distribution with d degrees of freedom (see Sections 7.2 and 7.3), and

.
I n the case of the gamma distribution model, the mHV

- DL is the x which satisfies 0.95 D **γ =** LCMRL Technical Report December 2010

Page 22 of 24 

### ( ) ( ) 

# ( **) ( ) ( )**  

# ( **)** Pr , 1 , C CD 

### Y L x x FL x x **αβ αβ γ ≥ =**  **− =** (20)

where 

### ( ) x **α** and 

### ( ) x **β** are defined in equation 19, and 

## ( )

, F **α β**

**⋅** is the gamma cumulative distribution function with parameters α and β

. For the data shown in Figure 9 -2, the mHV-DL and L c values are 0.026 ug/L and 0.016 ug/L, respectively. In th is case, which is typical, the mHV -DL is higher than L c , as one would intuitively expect. In the case of the TCE data used for Figure 9 -3, the reverse is true. The mHV -DL and L c values are 0.018 ug/L and 0.01 9 ug/L (both are converted from ng/L, as presented in Figure 9

- 3) , respectively.
Although this seems odd at first glance, in fact, the mHV -DL is in the domain of spiking concentrations while L c is in the measurement domain. Systematic low level bias in the measurements can cause this apparent anomaly. Even if a calibration does not show a statistical significant intercept term, because the LCMRL study involves much more data and because it includes preparation steps, it tends to show low level bias and nonlinearity not evident in the calibration data.

Let x c denote the spiking concentration at which the expected response equals L c . Then it is always true that x c < mHV -DL. In fact, x c is a very important quantity, since it is the natural
censor ing point to be used in data analysis for nondetects (measurements which are less than or equal to L c ). In the case of the TCE data used for Figure 9 -3, the intercept for the mean response model, at about 0.009 ug/L, is about half L c (0.019 ug/L), which i s a significant fraction of L c . Remember that L c is the estimated 95th percentile of the response distribution at zero spiking concentration. To get a response greater than 0.0 19 ug/L 95% of the time, we estimate that we must spike at 0.018 ug/L (the m HV -DL). Although it may seem odd, this situation can be explained by the positive intercept of the mean response model being “not small ” relative to L c .



## 11.0 Summary, Conclusions and Recommendations The proposed LCMRL procedure described in this report provides a statistically robust method of estimating the lowest concentration at which a laboratory can reliably achieve the MQO of 50%-150% recovery. Although the procedure is complex and computationally intensive, implementing it in user

- friendly software in th e public domain and freely available over the internet makes it easy for laboratories to use.
Furthermore, a robust method of estimating the LCMRL makes it possible, through Bayesian bootstrap resampling, to estimate by simulation the statistical distri bution of LCMRL , mHV

- DL and L c value s that would be generated by a randomly selected ‘experienced drinking water laboratory’ upon repeated execution of the LCMRL study design . This will enable the development of scientifically defensible MRL values for guidance and regulatory use.
Coupling the mHV

- DL procedure with the modeling used for the LCMRL computation improves the estimation of L c and the detection limit. Improvement of the estimation of L c will improve statistical inference regarding low

- level da ta sets. Improvement of the estimation of theLCMRL Technical Report December 2010

Page 23 of 24 detection limit will improve the accuracy of assessment of risk modificat ion due to regulatory actions.

The se development s have great potential value for regulators, the regulated community, risk
assessors and data analysts, but especially for laboratories. Use of these methods in the future should result in much better understanding and control of measurement data quality. In turn, this should lead ultimately to better data quality and decision-making.

The author has several recommendations for future work related to the LCMRL/Hubaux- Vos methodology. The first is the further characterization and potential improvement of the statistical methodologies proposed here, particularly viewing them as a combined procedure rather than isolated methodologies. Secondly is development of optimal design methodology (based on the statistical characteristics of the methodology) and development of more definitive recommendations for the design of LCMRL/DL studies.

The third recommendation is development or improvement of procedures for verification of MRL capability and ongoing monitoring of detection and quantitation capabilities in the context of a routine laboratory quality control program. The fourth recommendation is the development and release in the public domain of an R library (R Development Core Team, 2007) for t he LCMRL/Hubaux-Vos methodology, which will facilitate the use and further development of this methodology in the research communities of statistics, chemometrics, environmental analytical chemistry and forensic and clinical chemistry.



## 12.0 References Currie, L. A. (1968). “Limits for Qualitative Detection and Quantitative Determination.” Analytical Chemistry , Vol. 40, pp. 586-593. Dobson, A.J. (1983). In troduction to Statistical Modeling . Chapman and Hall, New York. Draper, N. and H. Smith (1981). Applied Regression Analysis: 2 nd Edition . John Wiley and Sons, New York.
Glaser, J.A., D.L. Foerst, G.D. McKee, S.A. Quane, and W.L. Budde (1981). “Trace Analy ses for Wastewaters.” Environmental Science and Technology . Vol. 15, pp. 1426-1435. Hoaglin, D.C., F. Mosteller, and J.W. Tukey (1983). Understanding Robust and Exploratory Data Analysis. John Wiley, New York.

Horn , P. S. (1988). “ A Biweight Predict ion Inte rval for Random Samples.” Journal of the American Statistical Association . Vol. 83, No. 401. (Mar., 1988), pp. 249-256. Hubaux, A. and G. Vos (1970). “Decision and Detection Limits for Linear Calibration Curves.” Analytical Chemistry . Vol. 42, No. 8, pp. 849-855.LCMRL Technical Report December 2010

Page 24 of 24 Jaynes, E. (1983). Papers on Probability, Statistics and Statistical Physics. A reprint collection

. R. D. Rosenkrantz (Ed.). Reidel. Dordrecht , Germany

. Kagan, A. M. ; Y.V. Linnik, and C. R. Rao (1973). Characterization Problems in Mathematical Stati stics . John Wiley. New York . 499 pp. Kelley , T. L. (1921). “ A New Measure of Dispersion.” Quarterly Publications of the American Statistical Association . Vol. 17, No. 134. pp. 743-749.
Lax, D. A. (1985). “Robust estimators of scale: Finite

- sample performance in long- tailed symmetric distributions ,” Journal of the American Statistical Association. Vol. 80, pp. 736-741. McCullagh, P. and J.A. Nelder (1989). Generalized Linear Models: 2 nd Edition . Chapman & Hall, London. Nelder, J.A. and R. Mead (1965). “A Sim plex Method for Function Minimization”, Computer Journal . Vol. 7, pp. 308-313. R Development Core Team (2007). R: A language and environment for statistical computing . R Foundation for Statistical Computing, Vienna, Austria. ISBN 3 -900051-07-0, URL http://www.R -project.org

. Randles, R.H. and D.A. Wolfe (1991). Introduction to the Theory of Nonparametric Statistics

. Krieger Publishing, Malabar, Florida. Rocke, D.M. and S. Lorenzato (1995). “A Two -Component Model for Measurement Error in Analytical Chemistry

. ” Technometrics . Vol. 37, No. 2, pp. 176-184. Rubin, D.B. (1981). “The Bayesian Bootstrap.” Annals of Statistics , 9, pp. 130-134. Thiel, H. and D.G. Fiebig (1984). Exploiting Continuity: Maximum Entropy Estimation of Continuous Distributions . Ballinger Publishing Co. Cambridge, MA. 246 pp. USEPA (1987). “National Primary Drinking Water Regulations- Synthetic Organic Chemicals; Monitoring for Unregulated Contaminants; Final Rule.” Federal Register . Vol. 52, No. 130. p. 25690, July 8, 1987.

USEPA (1999). “Revisions to the Unregulated Contaminant Monitoring Regulation for Public Water Systems; Final Rule. ” Federal Register . Vol. 64, No. 180. p. 50556, September 17, 1999.

USEPA (2004). Statistical Protocol for the Determin ation of the Single-Laboratory Lowest Concentration Minimum Reporting Level (LCMRL) and Validation of Laboratory Performance at or Below the Minimum Reporting Level (MRL). USEPA, Office of Ground Water and Drinking Water, Standards and Risk Management Division, Technical Support Center . EPA 815- R-05-006. Nov ember 2004.LCMRL Technical Report - Appendix A December 2010

Page A -1 

## Appendix A: Computation of MRL The computations carried out for establishing the Minimum Reporting Level (MRL) are programmed in the R language. The Matlab code for the LCMRL calculator was first translated into R and then amended to also produce estimates of the MRL estimate for all analytes appropriate for use in the LCMRL calculator. There are two reasons for programming the MRL in R; the first being that the recoding of t he Matlab code provided a validation check of the LCMRL code , and secondly it provided an interactive, non-compiled program that could be used internally by EPA to calculate MRL estimates and associated graphics.
The MRL is determined using a Bayesian bootstrap (BB) (Rubin, 1981) of the LCMRL estimator using the LCMRL study data from each of several experienced drinking water laboratories. The BB replicates that were generated from each laboratory’s data serve to estimate the distribution of estimated LCMRL values that each laboratory might generate on repeated performance of the LCMRL study . The distribution of pooled BB replicates, generated from LCMRL study data, can be used to approximate the distribution of estimated LCMRL values which might be generated from the p opulation of experienced drinking water laboratories.
The MRL is calculated in three steps whenever there are three or more laboratories providing data with valid LCMRLs or calculated LCMRLs that are below the lowest non -zero spiking level. In the first step, 200 BB LCMRL replicates are calculated for each laboratory data set. In the second step a predicted distribution of some unknown and yet to be observed laboratory is built from the population of replicate laboratory LCMRLs using a random effects model. In the third and last step the MRL is taken to be the upper 95% one-sided confidence interval on the 75 th percentile of the predicted distribution referred to as the 95-75 upper tolerance limit (95-75 UTL).
A description of the procedure is given in the following three steps.
Step 1: Bayesian Bootstrap of the Laboratory data When three or more laboratories have either a valid LCMRL or an estimated LCMRL that is reported to be less than the lowest reported non-zero spiking level a BB is carried out on each laboratory’s data to produce 200 replicate LCMRL estimates, l kr . Here k is the number of laboratories ( k **≥** 3), and r is the number of replicates (200)

, . The BB is implemented by selecting random Dirichlet weights from the uniform distribution on the n i – 1 dimensional simplex, where

n i is the number of measurements at the i th spiking level. These weights are associated with the original replicate measurements within each of the i spiking levels of the laboratory data. The LCM RL procedure is then rerun using the Dirichlet weights to compute a weighted LCMRL. This is done repeatedly to generate the sample of BB replicate LCMRL estimates. Using Efron’s bootstrap on a small data set, like the set of replicates at each spiking level, is not appropriate. Statistics computed using Efron’s bootstrap are essentially calculated as randomly weighted statistics where the set of possible weights for each observation is 

### { **}** 3 1 1 2 0, , , , , ,1 n nnn n **−**

 with the additional constraint that all o f the weights sum to 1. When n is small, the number of possible weights and therefore the number of values of any statistic that could be computed from the data given those weights, is small. Therefore our approach using the Bayesian Bootstrap with Diric hlet weights allows an extremely large number of weighted LCMRL values to be computed even from a relatively small data set.LCMRL Technical Report - Appendix A December 2010

Page A -2 The portions of the LCMRL calculation that use robust weights are modified to use a product weight of the form

*

* * 1 * , 1, , / , 1, , i D RW ij ij ij i

n ij ij ij i j www j n w w wj n

**= =**  **=**

**= =**  

## ∑ 

 where D ij 

### w is the Dirichlet weight and RW ij 

### w is either the robust Huber or Tukey bi square weight

. The weights are then normalized to produce the final weights w ij . Step 2: Estimating the Predictive Distri bution The predictive distribution of an unknown and unobserved laboratory from the random sample of laboratories reporting data is esti mated in three steps as follows

1. Power Transform of BB replicate Laboratory Data The BB replicate results are transfor med to approximate normality using the power transform (Box and Cox, 1964). The transformed replicate LCMRL values, , are cal culated as follows where **λ** is the transformation parameter

* , 1, , 1, ,200 if 0 ln( ), 1, , 1, ,200 if 0 kr kr kr l k Kr l l k Kr **λ**

**λ λ**  **λ **  **=... =... > ** 

**= ** 

**=... =... = **  **** where K is the number of laboratories. The transformation parameter λ is constrained to be non- negative.

2. The values of the predictive distribution are calculated according to the following step s:
a. Calculate Initial Resistant Location and Scale Estimates Location and scale parameters are estimated for the pooled set of transformed replicate values. These values are used as initial estimates in generating one

sided Tukey bi square weights to down-w eight large estimates i

* kr l n the pooled transformed replicate LCMRL values * kr l . The median of the is used for the initial resistant location estimate and the median absolute deviation from the median (MAD) of the * kr l is used as the initial resistant scale estimate.

b. Calculate One

- Sided Tukey Bi square Weights for the pooled LCMRL replicate values The one sided Tukey bi square weights are calculated using a tuning constant, c BW , of 6. The in itial starting values used are the median, , and MAD from “Step 2 a” above. If any of the weights are zero, the weights are recalculated with a tuning constant of 9. The weights are iteratively calculated at the t th stage accordin g to the following method * kr l

**( )** 0 BW 

### TLCMRL Technical Report - Appendix A December 2010

Page A -3 **( ) ( )**  **( )**
**( )**

**( )** 

# ( ) 

# ( ) **( )**

**( )**
**( )** 

# ( **)** 1 * 1 *

1 *

2

2 ( )

( ) ( ) , 1, , , 1, ,200, if 0, if
1 , if 1 , 1, , , 1, ,200 0, if 1 t kr BW

BW t l T
t kr BW c AD kr t kr BW

t t kr kr BW t kr
t kr

t BW t BW t BW kr kr kr kr kr k Kr l T u l T

u u w k Kr

u T wl w **−**

**− −** 

**− ** 

**= = >**  ****

**= **  ** ≤ **  ** − <**  ****

**= = =** 

** **  **≥ ** 

**=** 

## ∑ ∑  

 

u ntil a relative convergence of 1E

- 6 is reached for **( )** t BW 

### T

. After convergence the final one- sided robust weights are calculated as the normalized weights
( )

( ) BW t BW kr kr BW t kr kr w w w **=** 

## ∑

c. Calculate Robust means and Variances for each of the K laboratory replicate LCMRL sets Robust location and scale estimates are calculated by laboratory for the transformed replicate LCMRL values, * kr l

. These values are used in esti mating the predicted laboratories variance on the transformed scale. For the k th laboratory, an estimate of the mean, * k m , and variance , *2 k s , are calculated using the blended Huber and Tukey biweight methodology used in the LCRML routine except that a tuning constant of 6 is used instead of 9.

d. Estimate the variance of the p redicted d istribution Robust location and scale estimates are estimated by laboratory for the transformed replicate LCMRL values. The predicted variance, 2 pred 

### s , is then estimated as 

# ( ) **( ) ( )**  **( )** 

# ( ) **( )** 200 * 1 1

200 2 *2 1 1
2 200 * * 1 1 2 2 200
1 1
2 2 2 *

1
1 1/ w

b

b w K kr k k r

K kr k k r

K kr k k r

K kr k r

pred m w m

s w s

w mm s
w
s Ks s **= =**

**= =**

**= =**

**= = =** 

**=**

**−**

**=**

**−**

**=+ +** 

# ∑ ∑ **∑ ∑ ∑ ∑**  **∑ ∑**LCMRL Technical Report - Appendix A December 2010

Page A -4 e. Estimate the values of the predicted distribution The values of the predicted distribution are generated as follows. First the transformed values are generated as

* * **

, 

### * ( ) / , 1, , , 1, ,200 pred kr pred kr k k 

### l m s l m s k Kr **=+ − = =**

 
These values are then transformed back to the original scale as

*

,

, *

, exp(ln( ) / ), if >0 exp( ), 0 pred kr pred kr pred kr l l l if **λ λ λ**  ** ** 

**= ** 

**= **  ****

St ep 3 : Estimate the MRL as the 95

- 75 UTL of the Predict iv e Distribution The MRL is then estimated as the 95 -75 UTL from the values of the predicted distribution using a weighted version of the Guttman non-parametric procedure (Guttman, 1970) using the weights

kr w
.

References
Box, G. E. P. and D. R. Cox (1964). “An Analysis of Transformations,” Journal of the Royal Statistical Society , Series B . Vol 26, pp. 211-46.

Guttman , I. (1970). Statistical Tolerance Regions: Classical and Bayesian

. Hafner Press, Darien, CT. Rubin, D. B. (1981). “The Bayesian Bootstrap.” Annals of Statistics , 9, pp. 130-134.