Statistical methods for the blood beryllium lymphocyte proliferation test.

The blood beryllium lymphocyte proliferation test (BeLPT) is a modification of the standard lymphocyte proliferation test that is used to identify persons who may have chronic beryllium disease. A major problem in the interpretation of BeLPT test results is outlying data values among the replicate well counts (approximately 7%). A long-linear regression model is used to describe the expected well counts for each set of Be exposure conditions, and the variance of the well counts is proportional to the square of the expected count. Two outlier-resistant regression methods are used to estimate stimulation indices (SIs) and the coefficient of variation. The first approach uses least absolute values (LAV) on the log of the well counts as a method for estimation; the second approach uses a resistant regression version of maximum quasi-likelihood estimation. A major advantage of these resistant methods is that they make it unnecessary to identify and delete outliers. These two new methods for the statistical analysis of the BeLPT data and the current outlier rejection method are applied to 173 BeLPT assays. We strongly recommend the LAV method for routine analysis of the BeLPT. Outliers are important when trying to identify individuals with beryllium hypersensitivity, since these individuals typically have large positive SI values. A new method for identifying large Sls using combined data from the nonexposed group and the beryllium workers is proposed. The log(SI)s are described with a Gaussian distribution with location and scale parameters estimated using resistant methods. This approach is applied to the test data and results are compared with those obtained from the current method.


Introduction
Chronic beryllium disease (CBD), a disorder one of several criteria for diagnosis of the that mainly affects the lung, occurs in a disease (1). In vitro proliferation of small percentage of persons exposed to bronchoalveolar lavage (BAL) cells when beryllium dusts. Most investigators require exposed to beryllium is extremely sensitive evidence of beryllium hypersensitivity as to and specific for the diagnosis of CBD but is not suitable for screening since it is an invasive procedure (1). A noninvasive procedure based on the proliferative response of blood cells to beryllium has been developed and is referred to as the beryllium-specific lymphocyte proliferation test (BeLPT) (2). This modification of the standard lymphocyte-proliferation test is used to identify relatively rare individuals among worker cohorts who display delayed hypersensitivity reactions when exposed to beryllium metal. The BeLPT involves in vitro challenge of peripheral blood lymphocytes with salts of beryllium combined with assays for clonal proliferation of sensitized subsets of CD4 lymphocytes using tritiated thymidine uptake as a quantitative measure of blastogenesis. The test is conducted using 96-well microtiter plates; the amounts of tritiated thymidine incorporated by replicate wells containing lymphocytes challenged with beryllium is compared with uptake of radioactivity by replicate wells of nonchallenged lymphocytes to establish stimulation indices (SIs) as a measure of in vitro sensitivity to beryllium. A major problem in the interpretation of BeLPT test results is outlying data values ( 7%) among the replicate well counts.
The increasing use of beryllium in several new economic sectors emphasizes the need for medical surveillance in the workplace for CBD (3). In particular, beryllium has been used in the nuclear industry for a number of years. Kreiss et al. (4) have examined the epidemiology of CBD in a stratified sample of workers at a nuclear weapons plant and discuss the role of the BeLPT in beryllium disease surveillance in the nuclear industry. The U.S. Department of Energy (DOE) is operating a screening program for CBD that will eventually include approximately 15,000 current and former beryllium-exposed workers at 20 DOE sites. Each participating beryllium worker will have a BeLPT at an approved laboratory using a standard protocol developed by the Committee to Accredit Beryllium Sensitization Testing (CABST). The results of each assay will be evaluated and classified as normal, abnormal, or unsatisfactory.
A major concern that was not completely resolved by the CABST was how to deal with outliers that occur in the BeLPT data. The main purpose of this report is to propose a new statistical approach that can be used for analysis of a BeLPT assay that may contain multiple outlying well counts. Given their undue influence on the estimates of the SIs. a method for handling Environmental Health Perspectives * Vol 104, Supplement 5 * October 1996 outliers is needed. The current approach (as described in the July 1993 version of the CABST protocol ["Appendix"]) is based on an ad hoc outlier rejection method. As an alternative we propose using resistant estimation methods that are not sensitive to outliers. The BeLPT assay is described with a regression model that relates the expected well counts at each of the three beryllium concentrations to the control well counts for cells that are harvested after 5 and 7 days. Resistant fitting methods are used to estimate the SI for each of the six beryllium concentrations.
The main advantage of this approach is that estimates of the SIs are calculated without explicitly identifying and deleting the outlying well counts. A second question considered is the identification of beryllium-exposed workers who exhibit beryllium hypersensitivity. Most (over 90%) of the beryllium workers will have SIs similar to those of a control group with no known exposure to beryllium. However, even after the use of resistant estimation methods to minimize the effect of outlying well counts, the BeLPT for some beryllium workers will yield large SIs. In this case we want to identify the outliers (i.e., individuals with large SIs), since they represent beryllium workers who exhibit beryllium hypersensitivity.

Beryllium Lymphocyte Proliferation Test
A detailed description of lymphocyte culture methods, quality control measures, and examples of plate maps and printouts of raw data are included in the Appendix. Following is a brief description (Figure 1) of the protocol for the BeLPT culture assay as established by CABST and implemented by the BeLPT laboratory at Oak Ridge Institute for Science and Education (ORISE) as of July 1993. The details of this procedure and the equipment used vary at different laboratories that are performing the BeLPT.
First, a 15-ml blood sample is obtained from each patient and mononuclear cells are separated using density gradient centrifugation. Next, lymphocytes are cultured using standard methods at a final concentration of 2.5 x 105 cells per well in 96well, flat-bottom microtiter plates. For each BeLPT assay, 12 replicate control wells and four replicates for each experimental condition (i.e., 1, 10, and 100 pM of BeSO4, and mitogen-stimulated positive controls) are set up. Third, cells are incubated at 37°C for 5 and 7 days and a Table 1. Well counts for BeLPT assay (AC153 data shown pulse of tritiated thymidine is delivered prior to harvest. Cells are harvested on filter paper and counts are measured in a Packard Matrix 96 gas ionization counter (Packard Instrument Co., Downer's Grove, IL). Each filter is counted for 30 min and the results organized as shown in Table 1 for statistical analysis.

Statistical Methods
In this section we describe three methods of analyses for the BeLPT assay. The first method is the outlier rejection procedure proposed by CABST. A regression model is proposed to describe the two new methods, which are least absolute value (LAV) regression on the log counts and resistant maximum quasi-likelihood estimation. The regression model is motivated by the fact that the SI describes the relative increases in the proliferative response of berylliumstimulated cells to control cells. This leads to the log-linear representation of treatment effects. It is also apparent (Results) that the variability in the well counts increases approximately in proportion to the square of the expected count, so that the coefficient of variation (CV) is constant (i.e., the standard deviation is proportional to the mean). This implies that taking logs of the well counts leads to constant variance and additive effects, and the main parameters of interest are the log(SI)s. If there were no problem with outliers, standard least squares methods could be used on the transformed data. This approach was not considered since the occurrence of multiple outliers has been well established.  1   965  1173  828  862  Control  1  1474  7237  1021  976  Control  1  1500  1729  1672  1992  Be 1  2  1050  706  1434  687  Be 10  3  1551  1466  1661  2301  Be 100  4  3571  5780  4011  5229  Day 7  Control  5  9202  5253  3786  5212  Control  5  2310  2844  1915  3102  Control  5  2458  3936  3087  6588  Be 1  6  714  1135  6084  1097  Be 10   7  786  846  2757  652  Be 100  8  6037  8349  6852  10449  Day 5   PHA  9  82425  52954  52669  50487   Candida  10  35501  21623  21551 4 5 ConA (10 pM) 4 5 Ill. Harvest method (days 5  that used an ad hoc procedure for deleting outliers based on the value of the CV for each set of culture conditions. If the CV is greater than 0.3, the most extreme count is deleted. This procedure is continued until the CV is less than or equal to 0.3, provided no more than one-third of the well counts have been deleted. A patient's data are listed as "acceptable" if the resulting CV is less than 0.3 for both day 5 and day 7 control data, and for at least four of the six sets of beryllium-stimulated quadruplicates. If these conditions are not met, the BeLPT is called "unsatisfactory" due to high variability in the data. A BeLPT is also considered unsatisfactory if the positive-control response is too low (indicating lack of cell viability), or if the control counts are considered to be either too low (relative to background) or too high. We assume here that BeLPTs that are unsatisfactory for either of the latter two reasons are identified before further analysis using criteria that depend on laboratory experience.
The SIs for the stimulated cells are the ratios of the treatment means and the corresponding control means (after the outliers have been deleted), i.e.,

SI-mean(treated) mean(control)
The positive control wells are only counted for 10 min, so the SIs are multiplied by 3 to adjust for the counting-time difference. The results of applying this procedure to the BeLPT data in Table 1 are given in Table 2. These data are acceptable since both of the control CVs are less than 0.3, and all six beryllium-stimulated CVs are less than 0.3. The procedure used to determine if a BeLPT is abnormal is presented at the end of this section.

Rerssion Model for the BeLPT Data
Let Yjk denote the well count for the kth replicate of the jth set of culture conditions.
The expected count in each well is represented by a log-linear regression function: E(yjk) = Aj = exp(Xj13), [1] where j = 1. 0 and k = 12 for the controls and k = 1,2,3,4 for the berylliumstimulated cells and the positive controls.
In Equation 1, X1 is a row vector of indicator variables and / is the vector of regression parameters (below). We further assume that the variance of the well counts is proportional to the square of the In the absence of outliers, ordinary least expected count: squares on the transformed data would yield consistent estimates for the log(SI) parameters (5). The effect of outliers is Var(Yjk)= (OAj) [2] minimized by using LAV (or some other robust method) on the zjk. LAV regres-Equations 1 and 2 together are referred to sion-also known as LI norm, least as a generalized linear model with constant absolute deviations (LAD), and minimum coefficient of variation-as detailed by sum of absolute errors (MSAE)-is well McCullagh and Nelder (5). The distinct known to be resistant to outliers and is an values of the row vectors of covariates Xj, important particular case of a general class j = 1.10 are shown in Table 3.
of robust methods known as M-estimators With this parameterization, the first (6,7). In general, LAV regression requires three ,Bs represent the log of the SIs for the special computational resources to calculate three concentrations of BeSO4 on harvest parameter estimates (8). In this situation, day 5 and the next three Ps are the corre-however, it is only necessary to find the sponding estimates on day 7. The last two median of the log of the well counts for B s are the log(SI)s for the positive control each set of design conditions (say Z) and wells, and 07 and /8 represent the log of then subtract the control median for each the control well counts on days 5 and 7. harvest day from the beryllium-stimulated respectively. We have developed two out-medians. Frome et al. present details in lier-resistant approaches for estimating the Appendices A and C of a report for Oak SIs and the coefficient ofvariation, 4.
Ridge National Laboratory (9). A resistant estimate of the coefficient of variation can Se iond Method, Based on LAV then be obtained as Regression on Log(y) The first approach based on the regression = C X medianJ]Zmodel is to take the log of the counts since   Gaussian error model and for consistency with the usual least squares results, in which the estimated variance is multiplied by the correction factor n/(np) (10) and S-PLUS function MAD (11), which computes the median absolute deviation (MAD) estimate of the standard deviation. Alternative approaches to estimating 4 have been discussed in the context of LAV regression (7,12) and there is no consensus as to the best approach. In addition to the fact that 4 is of direct interest, it is also needed to obtain an estimate of the parameter covariance matrix (02(X'X)-l where ()2 = [2f(0)]2 is the asymptotic variance of the sample median (13). Following the approach of McGill et al. (14) we assume that the underlying error distribution is Gaussian in the center and use 6 = ,rI2~P to obtain an estimate of the standard deviation of the log of the stimulation indices. The appropriate diagonal term from (X'X)-l is 4/12, and consequently the estimated standard deviation of log(SI) is 1.25,L(O.58) = 0.72 4OL The results of applying this approach to the data in Table 1 are shown in Table 4.

Third Method, Based on Quasi-likelihood Estimation
In the second regression model approach, the analysis is done on the original scale and estimation is based on an iterative weighted least squares (IWLS) algorithm. The use of IWLS for generalized linear (15) and nonlinear (16) regression functions leads to maximum likelihood estimates when the dependent variable is in the regular exponential family. McCullagh (17) extended this result to quasi-likelihood (QL) estimation, which requires specification of the mean and variance function. Extension of IWLS to resistant/ robust regression has been described by Green (18) and Pregibon (19), and the computational approach described by Chambers and Hastie (20) is used here. Similar resistant regression methods have been applied to the analysis of drug concentration-time data encountered in human bioavailability studies (21). Consider the following weighted sum of squares, To adjust for the effect of outliers we introduce a second weight for each observation, -=f 1 lul<k W kl lul lul>k [5] [3] where Ai = exp(Xj,1) and Wjk a l/var(Yk) = l/2. The IWLS procedure starts with an initial estimate, say ,B°, and A7 in Equation 3 is replaced with the first order Taylor series exp(Xj13') + PjS' where P = XJ/, and the weights are evaluated at ,6 to obtain SWkYjk_(AO + pja) The unknown "correction vector," 60, is then calculated using weighted least squares, i.e. by solving [4] for 6'. The estimate of /B is then updated 13l = p + &0, and the procedure is repeated where u = (Yjk -3jk) I4)ik is the standardized residual using the current estimates of ,B and 4. This is known as an M-estimator with Huber's loss function. The "tuning constant," k, must be specified and we use k = 1.345, which leads to estimates with approximately 95% efficiency (19). Therefore, to obtain resistant quasi-likelihood estimates we multiply the elements of the diagonal matrix Win Equation 4 by the Huber weights discussed in Equation  5; again, details can be found in Appendix D in Frome et al. (9). Following the last iteration, an estimate of the coefficient of variation is obtained using a scaled MAD estimate of the standardized residuals Ujk= (Yjk -Yjk) jk 4 = 1.48 xmedian{lujIk} x n p Identification ofBeLPTs with Large SIs The CABST method for identifying an "abnormal" BeLPT is based on the  [6] where NA is the number in the nonexposed group. The mean and standard deviation of the MSIs are then used to calculate SI* = mean + 2 (standard deviations). A BeLPT for a beryllium worker is defined as abnormal if at least two SIs exceed SI*. Using the current value (SI* = 5.65), we conclude that the BeLPT in Table 2 is normal. The probability of obtaining a statistical false positive for this procedure is unknown.
We propose an alternative approach that establishes a reference database of BeLPTs based on BeLPTs from nonexposed workers and historical data from beryllium workers. The best way to establish this reference database is laboratory dependent and will not be discussed here. For illustrative purposes we will use the combined data from the beryllium workers and the nonexposed (control) group. The method is based on the assumptions that the estimates of the log(SI)s are approximately normally distributed and that almost all of the beryllium workers are not sensitized. Resistant methods are then used to counter the influence of outliers (i.e., the abnormal test results). The first step is to calculate an outlier-resistant estimate of location, ft, and spread, Sr, for each of the six log(SI) distributions in the reference database. In the results that follow we use j= median(1,1i, i =1.N), and i =MAD(13,1i = 1...,N), where N= 173 and j= 1, . . . 6. The second step is to convert the log(SI)s for each individual into standardized deviates ui.= - [7] S.
using the values of f1 and sL, from the reference database. The six standardized deviates for a BeLPT are compared to the zp quantile of the standard normal distribution. If at least two of these values exceed zp, the BeLPT is called abnormal. If the estimated log(SI)s are independent, then the binomial distribution can be used to calculate an approximate probability of at least k out of six "large" SIs for a given value of zp.
The probability of at least one large SI is 1 -p6 = 0.141 (forp = 0.975). The probability of at least two is 1 -[p6 + 6(1 -p)p5] = 0.009 (for p = 0.975). In fact, the log(SI)s are positively correlated, so this probability should be a lower bound on the chance of finding a false positive BeLPT.

Results
The regression model and the estimation methods were obtained through analytic reasoning and limited experience with a few data sets. To evaluate the utility of our two new methods, we applied them to all available BeLPT assay results obtained at the ORISE BeLPT laboratory as of July 1993. The outlier rejection method in use at ORISE at that time also was applied to each BeLPT assay.
As a preliminary step we provide an analysis of the 12 replicate control wells on days 5 and 7 for each of the 173 BeLPTs. Estimates of scale and location are computed to verify the assumed form of the mean-variance relation. We then describe the distribution of estimates of the log(SI)s for each beryllium concentration. Under the null assumption that if an individual is not sensitized to beryllium his/her SIs should be one, and the estimates of the log(SI)s will be approximately normally distributed with a zero mean and the covariance matrix indicated in "Methods." The true SI for an individual for any given beryllium concentration is, of course, unknown. In a population of nonsensitized individuals, the true log(SI)s may differ from zero. Consequently, the distribution of estimates of the log(SI)s presented in this 'section reflect sampling variation, possible differences in responsiveness among individuals who are not beryllium sensitized, and the presence of beryllium-sensitized workers. As a matter of convenience, we may refer to the distribution of estimates of the log(SI)s in this section as a distribution of log(SI)s.

Description of the Data
A total of 173 BeLPTs are used in this evaluation. There are 133 from a group of 120 workers exposed to beryllium; the remaining 40 are from persons who have no known exposure to beryllium. The discrepancy between the number of test results and the number of beryllium-exposed workers is accounted for by the fact that a second BeLPT was carried out on 13 workers. Ideally, there should be 56 observations (well counts) for each assay, but in some cases, well counts are missing due to insufficient number of cells or technical errors ("Appendix"). When an assay is incomplete, parameters are estimated (if possible), based on the reduced data set. Comparison ofMoment and Resistant Estimates ofthe Coefficient of Variation for Control Wells An important assumption is that the standard deviation of the well counts is proportional to the mean as implied by Equation 2. Each of the 173 assays contains 12 replicate control wells on both days 5 and 7. To verify this assumption, location and scale estimates for the control wells for each assay on days 5 and 7 were calculated. Figure 2A shows the relationship between the moment estimator of location (y, the sample mean) and the moment estimator of scale (s, the sample standard deviation) for the day 5 control wells.
The solid line is the least squares fit s= 0.448y, and the slope (0.448) is an estimate of the coefficient of variation for day  Figure 2B reflects the use of resistant methods. Figure 3 shows the relationship between the resistant estimates of location and scale for the day 5 control wells ( Figure 3A) and the day 7 control wells ( Figure 3B) on a log-log scale. The solid line for day 5 in Figure 3 (A) is log& = log(O.34) + log(]), and the solid line for day 7 ( Figure 3B) is logd = log(0.36) + log(f). Note that if& is proportional toy (i.e., constant coefficient of variation) then the log-log plot should be linear with a slope of one. The slope of the least squares regression of log a on log(y) for day 5 is 1.04 (standard error = 0.04) and for day 7 the slope is 0.96 (standard error = 0.03). Since neither estimate is significantly different from 1, this supports the regression model assumption of constant coefficient of variation. The main difference in the day 5 and day 7 results is that the day 7 results are shifted to the right since the control-well counts are generally higher on day 7 than those on day 5. The median of the ys on day 5 is 1247 compared with 1840 for day 7. These results are consistent with the laboratory observation that day 7 results are generally higher and show greater variability than well counts on day 5.

Summary ofResults for The Methods
The three methods of analysis were applied to the data described at the beginning of this section. For each method, two graphical displays summarized the results. Only the results for the LAV method are presented here since the plots for the other two methods were very similar in appearance and are available elsewhere (9).
The first graphical display ( Figure 4) is a series of 12 boxplots (11,14,22) placed side by side for the log(SI)s; the horizontal axis on the bottom shows the untransformed SIs. The ends of the box correspond to the 25th and 75th percentile so that 50% of the log(SI)s are contained in the box for each group. The vertical dotted lines are drawn to the nearest value not beyond a standard span-1.5 x (interquartile range)-from the quartiles. The outlying values are shown individually for each group of data. There are two boxplots for each beryllium concentration on days 5 and 7. The first one in each pair is labeled -6 -4 -2 "BW" for beryllium workers, and the second one is labeled "NE" for not exposed. Consequently, each pair of boxplots provides a comparison of the distribution of the SIs for the beryllium group and the nonexposed group for each of the six culture conditions. Consider, for example, the first two boxplots in Figure 4 which are for beryllium concentration 1 on day 5 (BW-D5Bel and NE-D5BeI) for the LAV estimates. Both log(SI) distributions are centered near zero and the nonexposed group is a little more spread out in the center. The beryllium-workers group shows nine outlying values in the positive direction and one in the negative direction. The notches (which represent confidence limits for the sample median) in the boxplots overlap, indicating that the difference in the location of the two distributions is not significant at a rough 5% level. The broken vertical line corresponds to log(SI) equal to zero and passes through both notches, indicating that both distributions are centered near zero on the log scale. Each of the boxplots in Figure 4 is centered near zero and is spread out evenly in both directions.
As the beryllium concentration increases on each day the variability (as indicated by the length of the boxes) increases. On day 7 the SIs for the beryllium workers are generally Log (SI)  smaller that those for the nonexposed group, and the median log(SI)s are less than zero except for the NE-D7Be1OO group. Results on day 7 are more variable than those on day 5 for each of the three beryllium concentrations. We have no definitive experimental data to explain the larger variability on day 7. However, considering that we are conducting short-term cultures of lymphocytes without replenishing of medium, it is not surprising to observe greater variability among wells with increasing time in culture since one would expect to see some depletion of nutrients, accumulation of metabolic byproducts, or scenescence of lymphocyte over time. In beryllium-sensitive persons, increasing divergence of counts between replicate wells would be anticipated to result from clonal expansion of sensitized CD4 subsets which would be expected to become more pronounced with increasing numbers of cell replications. The second graphical display ( Figure 5) shows a normal (Gaussian) probability plot for the combined BW and NE SIs for each of the three beryllium concentrations on days 5 and 7 (23). In each of the six plots, the data (ordered values of the log(SI)s) are shown on the vertical axis on the left, and the quantiles of the standard normal distribution are shown on the horizontal scale. Statistical theory indicates that estimates of the log(SI)s should be approximately normally distributed, and the large sample standard deviation should be about 0.28 if the coefficient of variation is 0.4 (see "Statistical Methods"). If the relation between the empirical and theoretical quantiles is linear, this indicates that the distribution is Gaussian. In each plot we have included the median (labeled M) and a resistant estimate of the standard deviation (labeled S) for the log(SI)s. The solid line in each plot shows the relation expected if the log (SI) values are from a normal distribution with median M (which determines the intercept) and standard deviation S (which determines the slope).
(The values of M and S are also shown in Table 5). Resistant methods were used to estimate the location and scale parameters for the combined data from the BW and NE groups. This reflects the assumption that most beryllium workers do not show an abnormal response, i.e., they look like the nonexposed group. For example, consider the plot for day 5 Be 1 in Figure 5. The log(SI)s appear to be approximately normal in the center. but several values are larger than expected (these are the points above the line). These outliers are SIs that indicate hypersensitivity to beryllium.
The results in Figure 5 indicate that the log(SI)s are approximately normally distributed. The center of each log(SI) distribution is greater that zero for each beryllium concentration on day 5 and is less than zero for each beryllium concentration on day 7. The untransformed SI units are shown on the vertical scale on the right side of each plot. The estimated standard deviations increase with beryllium concentration on each day, and are larger on day 7 than on day 5. The estimates of location (,u) and scale (s) for each method are summarized in Table 5. For each beryllium concentration the estimates from the three methods are in very close agreement. The boxplots and normal probability plots for the current method and the QL method are not shown here since they are almost identical to Figures 4 and 5 and are available elsewhere (9).

Comparison of the Current Method and LAV Method
A direct comparison of estimates of the log(SI)s obtained using the current method and the LAV method for each beryllium concentration is given in Figure  6. The slope of the line in Figure 6 indicates exact agreement between the two methods. To further compare the current method and the LAV approach we use the average difference (avedif) of the log(SI)s. For each BeLPT the LAV log(SI)s are subtracted from the corresponding log(SI)s based on the current method. The result is  For example, for AC 147 (Table 6), the current method day 5 BelOO SI is exp(1.41) = 4.10 and exp(1.45) = 4.26 for the LAV procedure. The log percent (L%) difference is 100 x log(4.10/4.26) = 100 x (1.41 -1.45) = -4L% where L% stands for the logarithmic percent (24). The SI for the current method is 96% of the LAV SI, i.e., about 4% smaller. The average difference between the current method and LAV method for AC117 is -1.7L%. Table 6 compares AC147s log(SI)s for all three methods. Figure 7 (left panel) shows a boxplot of avedifl2 (as defined above) for the 173 BeLPTs. The average difference is between Table 6. Comparison of log (SI)s for patient AC147.a -3L% and 9L% 50% of the time and the mean of the average difference is 3L%. A large, positive value of avedifl2 indicates that the current method log(SI)s for a BeLPT are greater than the LAV log(SI)s. The results in Figures 6 and 7 and Table 5 show that the estimates of the log(SI)s for the LAV and current method are in very close agreement.

Comparison ofResistant Quasi-Likelihood Method and LAV Method
The plot ( not shown) of the QL and LAV log(SI)s was almost identical to Figure 6. The close agreement between the two methods is further demonstrated by the average difference avedif32 = mean[100 * (lBQL -'L0 j = 1.6].  where a large positive value indicates that the QL log(SI)s for a BeLPT are greater than the LAV log(SI)s. Figure 7 (right panel) shows a boxplot of these values. The average difference is between -6L% and 3L% 50% of the time, and the mean of the average difference is -0.6L%.

Identification ofBeLPTs with UV SIs
The first step in the alternative method is to convert each log(SI) into a standardized deviate (Equation 7) using the values of Uj and ij given in Table 5. These standardized deviates are compared with the quantiles of the standard normal distribution, i.e. Pr[u < zp] = p. If at least two of these u ,s exceed zp that BeLPT is called abnormal. The results (i.e. the u11s for the abnormal BeLPTs ) of applying this procedure to the 173 assays using zO975 = 1.96 are shown in Table 7. In particular. the standardized deviates for the LAV log(SI)s for AC147 (Table 6) are given in row 2.
An alternative to calculating standardized deviates for each log(SI) is to calculate critical value for each SI SI; = exp(,t7i + zp j), j =1.6.
The critical values obtained using fj and for the LAV method in Table 5 are shown in Table 8. If a patient's SI exceeds the corresponding critical SI, that SI is considered large. Thus, if any two SIs exceed the corresponding critical SI, that patient's data are considered abnormal. For example, the LAV SIs from Table 6 (patient AC 147) are 1.27, 6.36, 4.26, 0.99, 6.11, and 6.11.  lDs marked with t were identified as abnormal by the CABST method (Section 3.5) and those marked with u were unsatisfactory. One BeLPT (AC149) called abnormal by the CABST method does not appear in this table. Since the day 5 Be 10 SI (6.36) and the day 7 Be 10 SI (6.11) exceed their respective critical SIs, these data are deemed abnormal. Regardless of the formulation used (calculating the standardized deviates or comparing the SIs to critical SIs) the conclusion is identical. The CABST procedure currently used at ORISE to identify workers with large SIs is based on the distribution of the maximum SI for each individual in the nonexposed group (Equation 6). A BeLPT for a beryllium worker is defined as abnormal if at least two SIs exceed SI* (currently equal to 5.65). This leads to the identification of patients AC147, AC149, AC235, and AC236 as abnormal. Using the outlier rejection method, six BeLPTs that had two or more SIs greater than 5.65 were found to be unsatisfactory based on the values of the within-group CVs. They are patients AC161, AC171, AC174, AC196, AC208, and AC225. Three of the abnormal BeLPTs are listed in Table 7 (only AC149 is missing). Table 7 also lists five BeLPTs as abnormal that were not identified by the current method as either abnormal or unsatisfactory.
In situations in which there may be excess variability, the CV can be used to evaluate the quality of the BeLPT. For the LAV approach L is a resistant estimate of the "within-group" standard deviation of the log well counts. Since L is not inflated by a few outliers (that could be caused by measurement error), it may be reflecting some intrinsic biological variability associated with the lymphocyte proliferation response in certain cell donors.

Conclusions
Three approaches to the analysis of the BeLPT have been described. The first method is the outlier rejection procedure (in use at ORISE in July 1993), and two new methods (LAV and QL) are based on resistant regression techniques. Each method was applied to a database of 173 BeLPTs (133 from beryllium workers and 40 from individuals with no beryllium exposure). Graphical and numerical summaries show that the three methods are generally in very close agreement. Both of the new methods are highly resistant to outliers (in the well counts), have wellknown statistical properties, and provide a "pooled" estimate of the coefficient of variation (0) for each BeLPT. The QL  Figure 9. Normal probability plot of log ('L)method requires an iterative algorithm and does not appear to offer any practical advantage over LAV. The LAV method is also easy to understand and compute and is recommended for routine analysis of the BeLPT.
Estimates of the log(SI)s are approximately normally distributed. The log(SI) distributions are centered near zero for each of the three concentrations of BeSO4 on harvest days 5 and 7. The variability is greater on day 7 than on day 5, and increases with concentration on each day. Resistant estimates of the location and scale parameters for each of the six log(SI) distributions are used to define large SIs, which are used to identify abnormal BeLPTs. Results of this preliminary approach to identify abnormal BeLPTs were compared with results obtained using the current method, and the discrepancies between the two methods suggest that a more detailed evaluation of the procedure is needed.
In a subsequent report further consideration will be given to the use of the LAV approach to address the following questions: a) How should "abnormal" BeLPTs be identified? b) Should a BeLPT be considered unsatisfactory as the result of high variability? c) How should the resistant estimate of the coefficient of variation (OL) be used in the BeLPT analysis?
The methods developed will be applied to a much larger database of BeLPTs obtained from the ORISE BeLPT laboratory and at least one additional laboratory that is currently using this assay to identify persons who may have CBD. The data set used in this report and an electronic version of ORNL-6818 (9)  The ORISE protocol for performing lymphocyte proliferation assays essentially adheres to the recommendations of the expert panel (i.e., CABST) convened jointly by the U.S. DOE Office of Health and the Beryllium Industry Scientific Advisory Committee (BISAC) at a meeting held in Washington, DC, on February 3-4, 1992. We collect approximately 30 ml of venous whole blood in sterile vacutainers containing sodium heparin for each assay (Figure 1, text). Tubes are inverted to mix blood with the anticoagulant and transported to the laboratory for processing. Cells are maintained at room temperature overnight. Within 24 hr after blood collection. mononuclear cells are separated using Ficoll-hypaque density gradient centrifugation, carried through three sequential washes and counted in triplicate on an automated cell counter. Lymphocytes are cultured in RPMI 1640 culture medium (GIBCO Grand Island, NY) buffered with Hepes salts, and supplemented with 2 mM /-glutamine, 100 units/ml penicillin, and 100 pg/ml streptomycin. Pooled human serum is added at a final concentration of 10%. We are using 96-well, flat-bottom microtiter plates and a final cell concentration of 2.5 x105 cells per well contained in 0.2 ml volume of medium.
Beryllium sulfate ([BeSO4], Brush and Wellman, Elmore, OH, 99.9% purity) in concentrations of 1, 10, and 100 pM is used to evaluate donor lymphocyte hypersensitivity to Be metals. As positive controls we use concanavalin-A (10 pg/ml) and phytohemagglutinin (30 pg/ml). For each set of exposures, quadruplicate wells are evaluated to obtain estimates of lymphocyte proliferation response. Unstimulated control wells are run in replicates of 12 because other laboratories have observed considerable variability in rates of tritiated thymidine incorporated in the control series, and extra replicates are needed to achieve the required levels of statistical confidence. All cells are incubated at 37 + 0.5°C in an atmosphere of 5% CO2 in air. Cells assayed for response to Be are harvested at 5 and 7 days with a terminal 6-8 hr pulse of 1.0 pCi of tritiated thymidine (specific activity 6.7 mCi/mM). We are using a Packard 96-well cell harvester (Packard Instruments, Downer's Grove, IL) that deposits lymphocytes from each individual well on a standard glass filter paper; the lymphocytes then can be counted intact on the Packard Matrix 96 gas ionization counter, or punched for assay using a liquid scintillation counter. The Matrix 96 unit is less efficient in detecting beta decays than scintillation counters but has the great advantage of simultaneously detecting beta radiation emissions from all 96 wells. Statistical accuracy can be achieved quite readily by increasing counting time using this instrument.

Quality Control
Excess variability in counts between replicate wells within a treatment, i.e., outliers could result from technical errors in initiating the tests, or possibly from intrinsic biological variables associated with the characteristics of lymphocyte proliferation response in certain cell donors. Sources of technical error might include mistakes in pipetting, such as failures to add appropriate numbers of cells to individual wells, lack of addition or double addition of tritiated thymidine to specific wells, or improper washing of filters resulting in residual counts of unincorporated thymidine, or smearing of radiolabel across the filter paper.
Stringent methods for quality control are used routinely to guard against inadvertent technical errors. To minimize the risk of pipetting errors, all media and other test reagents are delivered to complete rows or columns of the test plate using electronic micropipetters that deliver up to 8 or 12 aliquots simultaneously. Thus, it is not likely that the operator could "loose her place" in adding reagents. Cells are harvested onto the surface of filter paper using the 96-well harvester, which simultaneously aspirates the cellular contents from each well. To ensure complete washings of culture plates, a wash volume of approximately 10 times that recommended by the manufacturer is used. For all tests, we routinely leave all wells in rows A and H empty as a quality control measure to allow evaluation of background counts on both the top and bottom of the filter paper. Erratic or high counts in these empty wells would signal incomplete washing of plates Environmental Health Perspectives -Vol 104, Supplement 5 * October 1996 or "smearing" or radioactivity from one well to another. Filter papers are counted intact on the Matrix 96 gas ionization counter, which simultaneously records counts and countsper-minute with attendant errors for each well. Because the matrix counter is a gas ionization unit, only those beta decays that are emitted at right angles to the surface filter pad are detected and recorded. Thus, the sensitivity of the instrument in detecting counts is considerably less than that of a liquid scintillation counter (about 20% of emissions are detected using the gas ionization unit). For this reason, all plates are counted for longer periods to accumulate enough counts for statistical accuracy. Routinely, all plates containing control wells and wells challenged by beryllium salts are counted for 30 min, whereas mitogen-stimulated positive controls are counted for 10 min each.
plate. Cells from patient 1 are pipetted into columns 1 to 4; cells from patient 2, into columns 5 to 8; and cells from patient 3, into columns 9 to 12. Rows A and H are left blank to monitor background counts in the culture system. Rows B, C, and D are replicate sets of control wells, whereas rows E, F, and G contain beryllium concentrations of 1, 10, and 100 pM, respectively. The lower half of the figure demonstrates the platemap for initiating cultures with phytohemagglutinin or ConA. An example of a typical printout of data from three different individuals is shown in Figure A2. The test is a 5-day plate, counted for 30 min. Data are shown as total counts. Patient 1 displays a pronounced response to all three levels for beryllium-challenged wells, whereas patient 3 demonstrates higher levels of counts in control wells but also demonstrates no response to beryllium. Direct comparisons of data among the three persons can be readily made from a single printout sheet. This approach allows comparisons of counts within replicate treatments for lymphocytes from the same donor, as well as comparisons of inter-individual variability in counts between different subjects.