UDC 613.632.4+614.72]:613.155.3-07:519.24

V. G. Kovyazin

Particular Features of Experimental Design and Multivariate Statistical Analysis in Determining Threshold and Harmless Levels of Exposure to Harmful Chemical Substances

S. V. Kurashov Kazan Medical Institute

One way to improve the effectiveness of scientific research on the establishment of hygienic standards for harmful chemical substances in occupational environments and in the biosphere may be the widespread use of modern computing technology and mathematical methods to determine threshold and harmless exposure levels.

In our view, a key aspect of this problem is the standardization of methodological approaches and the development of appropriate computational algorithms based on them. The use of even the most advanced computers in combination with traditional methods clearly cannot be regarded as a qualitatively new methodological approach.

A number of authors have repeatedly pointed out shortcomings in the existing design of sanitary-toxicological experiments [5, 9–11]. The main emphasis has been on the lack of a probabilistic assessment of threshold and no-effect doses (concentrations), and consequently their lack of accuracy, as well as the absence of clear diagnostic criteria for the state of the organism as a whole at different levels of toxic exposure. The question of reducing the duration and labor requirements of chronic experiments remains particularly pressing.

Despite their specific features, these shortcomings are typical of situations in which the levels of exposure applied to a system and its response are known, but information about processes within the system is limited (the “black box” principle). To solve problems of this kind, methods of mathematical experimental design [1, 4], among others, have been developed in recent decades and applied in various branches of science. The essence of these methods is as follows. It is assumed that the behavior of the system at different levels of exposure to the factors can be

described by a specific mathematical model. Using a predetermined design and a minimum number of experiments, mathematical statistical methods are used to establish the correspondence between the relationship actually observed and the theoretical model and to determine its specific parameters (the coefficients of the equation). The mathematical model makes it possible to predict by calculation both the maximum possible response level (the largest value of the function) and the relationship between the exposure levels of the factors under study at specified response values.

The methodological approach we propose is based on the existing framework for sanitary standard setting. It involves using, in the design process, a mathematical model of the dose–time–effect relationship, together with a comprehensive assessment of health indicators by multivariate statistical methods. The latter evidently brings us closer to assessing toxic effects at the level of the whole organism, while the mathematical model makes it possible to approximate the experimental data correctly and obtain a probabilistic estimate of the threshold for toxic effects and of no-effect levels. Essential requirements are the measurement of the entire set of indicators at every stage of the experiment and an unambiguous interpretation of the concept of a threshold dose. Hereafter, by the threshold for toxic effects we mean the minimum doses (concentrations) that cause statistically significant deviations in the indicators studied in experimental animals compared with controls. Given the controversial nature of “true homeostatic adaptation” and, above all, the absence of objective criteria for it, we regard improvements in the health of test animals as one of the phases of the toxic process [8] and take them into account, as we do deviations of indicators toward deterioration.

Let the doses (concentrations) tested in the experiment be denoted by \(X_1\), and the time at which the indicators of toxic effects are recorded by \(X_2\). The dose–time–effect relationship can then be represented in general form as a function:

\[ Y=f(X_1,X_2), \tag{1} \]

where \(Y\) is a variable dependent on \(X_1\) and \(X_2\) (the “response” of the mathematical model). \(Y\) is intended to reflect the health of the animals in the experimental groups at the level of the whole organism. As a first approximation, it seems logical to take \(Y\), by analogy with probit analysis, as the proportion (in %) or probit of “sick” animals, i.e., animals whose indicators fall outside the normal range. We previously described a probabilistic method for estimating the limits of the normal range for identification based on a single indicator [7]; in principle, there are no restrictions on applying it to several indicators. Preliminary calculations have shown, however, that such identification yields a coarser estimate: doses that differ significantly from the control when compared by conventional methods are classified as certainly having no effect. Such an “upheaval” in sanitary standard setting would be undesirable, since it would require revision of all calculation formulas and principles for extrapolating from \(\mathrm{Lim}_{ch}\) to the maximum permissible concentration. In view of the foregoing, it would apparently be appropriate to use as \(Y\) a statistical criterion—a multivariate analogue of Student’s \(t\)-test or Fisher’s \(F\)-test—that would make it possible to assess the significance of the control–experimental difference simultaneously across the entire set of indicators studied. One such measure, Hotelling’s \(T^2\) statistic, has been described [2, 3]:

\[ T^2=\frac{n_1n_2}{n_1+n_2} (\overline{X}_i-\overline{X}_j)'S^{-1} (\overline{X}_i-\overline{X}_j), \tag{2} \]

where \(\overline{X}_i\) is the vector of mean indicator values in the experimental group; \(\overline{X}_j\) is the corresponding vector for the control group; and \(S^{-1}\) is the inverse of the pooled covariance matrix \(S\):

\[ S=\frac{1}{n_1+n_2-2} \left[x_1'x_1+x_2'x_2\right], \tag{3} \]

where \(x_1\) is the matrix of observations in the experimental group (expressed as deviations from the mean indicator values); \(x_1'\) is the transpose of \(x_1\); \(x_2\) and \(x_2'\) are the corresponding matrices for the control group; and \(n_1\) and \(n_2\) are the numbers of observations. Despite the somewhat cumbersome appearance of formulas (2) and (3), expressed in the notation of matrix algebra, values of the \(T^2\) statistic can be calculated fairly quickly on an Elektronika D3-28 desktop microcomputer. Algorithms for calculating \(T^2\) and FORTRAN programs are described in the literature [3]. We have tested these programs and adapted them to this microcomputer. It should be noted that the number of indicators studied must always be less than \(n_1\) and \(n_2\); otherwise, the program will fail because of the characteristics of matrix computations.

Having defined \(Y\), \(X_1\), and \(X_2\) in (1), we proceed to establish the mathematical model. For this purpose, we turn to the computational algorithms and terminology used in the practice of mathematical experimental design [1, 4], and in particular to the algebraic polynomial model, which has been developed in detail in the applications described there. For our case of two independent variables \(X_1\) and \(X_2\), and taking into account the available data on the possible form of the dose–time–effect relationship in a chronic experiment [8], a sufficiently close approximation to the actual relationship can be

Calculation table for determining the polynomial coefficients

was obtained using an incomplete third-degree polynomial, corresponding to a two-factor \(3^2\) design.

In our case, unlike the problems described in [1, 4], a model with parameters that allow prediction of the maximum value of \(Y\) is not essential. Our aim is to determine the threshold dose, that is, the value of \(X_1\) corresponding to the maximum \(X_2\) and to the value of \(Y\) corresponding to the critical \(T^2\)—its minimum value at which the control–experimental difference is significant at the specified level of confidence in the statistical estimates.

To illustrate the foregoing, we present our data supporting the maximum permissible concentration of ethylene chlorohydrin in workplace air [6]. The table presents the initial data for calculating the coefficients of an incomplete third-degree polynomial. The values of \(X_1\) and \(X_2\) are presented in the table (columns 1 and 2) in coded form: −1, minimum; 0, intermediate; and +1, maximum. For \(X_1\), these correspond to 0.1, 1, and 10 mg/m³; for \(X_2\), to 1, 2, 3, and 4 months. (Four levels were used for \(X_2\), although three are sufficient under a \(3^2\) design: 1, 2.5 [code 0], and 4 months; these intervals, however, do not correspond to those normally used in sanitary-toxicological experiments.) The data in columns 3–7 were obtained by arithmetic operations according to the expressions in the table header. The data in column 8 were calculated using formulas (2) and (3) for two measures: body weight and the value of the CNS summation-threshold index (the other measures were investigated periodically or at the end of the experiment). Coding facilitates calculation of the polynomial coefficients and provides a convenient graphical response surface irrespective of the units of \(X_1\) and \(X_2\). Coding and conversion from coded to actual values are performed using simple formulas [1, 4].

The polynomial coefficients were calculated from the table data by the method of least squares, and our dose–time–effect model takes the form:

\[ Y=4+5.3X_1+1.29X_2+3.71X_1^2-1.27X_2^2 +3.62X_1X_2+1.78X_1^2X_2. \tag{4} \]

The polynomial term \(b_{122}X_1X_2^2\) was not included in the final model because \(b_{122}=0\).

The statistical measures of model adequacy are: multiple coefficient of determination \(R=0.947\), standard root-mean-square error of the response \(Z=2.06\), and \(F\)-ratio \(F=14.99\). The calculated values of \(Y\) are given in column 9 of the table.

It should be emphasized that model (4) is valid only for this example and must be justified anew each time a particular set of experimental data is analyzed. Moreover, where advanced computing facilities are available, it is advisable to calculate several models and compare their adequacy with that of the polynomial model using the statistical criteria mentioned above, primarily the \(Z\) criterion.

A mathematical model of the dose (concentration in our example)–time–effect relationship makes it possible to present the experimental data as a nomogram, a geometric representation of the response surface. The figure shows such a nomogram for our example. To construct it, approximate values were calculated on a microcomputer for \(Y=1\ldots18\) and \(X_2=-1\ldots+1\) in increments of 0.1

Dependence of the significance of the control–experimental difference on concentration \(X_1\) and time \(X_2\) in a chronic experiment. Geometric representation of the response surface (Hotelling’s \(T^2\) criterion). Two-factor model of an incomplete third-degree polynomial (see text for explanation).

The abscissa shows experimental duration (months); the ordinate shows substance concentration (mg/m³, logarithmic scale). Curves 1–7 are isoeffect lines for values of \(X_1\) and \(X_2\) corresponding to the respective \(T^2\) values (1–7). The shaded region represents values of \(X_1\) and \(X_2\) for which the control–experimental difference is statistically significant. The intersections of lines \(A\) and \(B\) with the ordinate \(A\) correspond to \(\mathrm{Lim}_{\mathrm{ch}}\), estimated at confidence probabilities \(P=0.05\) and 0.01.

…of \(X_1\) (by bisecting the interval). Points with coordinates \(X_1, X_2\) corresponding to equal values of \(Y\) were joined by smooth isolines. The resulting nomogram resembles a topographic map and provides a clear and accurate representation of the complex “response” surface. Of particular interest are the minimally effective concentrations that produce statistically significant differences from the control at the end of the chronic experiment. These are points \(A\), at the intersection of \(X_2=+1\) with the isoline \(Y=7.605\) (the critical value of \(T^2\) at \(P=0.05\), taken from the table of values for Fisher’s \(F\) test), and \(B\), at the intersection with the isoline \(Y=12.94\) (\(p=0.01\)). The numerical values of \(X_1\) at points \(A\) and \(B\) are 0.334 and 0.699, respectively. Taking the antilogarithms of these values (the numerical concentration values were logarithmically transformed before coding) gives threshold concentrations at the respective confidence levels of 2.1 mg/m³ (\(P=0.05\)) and 5.0 mg/m³ (\(P=0.01\)).

The algorithm described is to some extent universal and can, with minor modifications, be used in other studies of the dose–time–effect relationship: determining threshold doses in acute and subacute experiments; calculating the median effective time to death when determining \(\mathrm{DL}_{50}\), with probits as the “response”; assessing the combined action of two substances (in which case \(X_1\) and \(X_2\) are dose levels and \(Y\) is the recorded effect); optimizing the dose and frequency of drug administration in pharmacology and experimental therapy; and so on. A particular issue is the choice of dose (concentration) levels in a chronic experiment. In our opinion, it is advisable to recommend as the maximum level (code \(+1\)) the threshold dose (concentration) determined using the above algorithm and identical endpoints in a 1-month subacute experiment, and as the middle level (code 0) the calculated tentative safe exposure level (ОБУВ). The minimum level (code \(-1\)) can be determined as follows:

\[ \lg X_{\min}=\lg X_{\text{ОБУВ}}- \left(\lg X_{\max}-\lg X_{\text{ОБУВ}}\right). \tag{5} \]

A separate issue requiring special consideration is the choice of health-status indicators, since their number is limited by the size of the animal groups (see above) and, for example, cannot exceed 9 when \(n=10\). Several models can, however, be calculated: for integrated indicators and specific tests, and even for individual tests (in which case Student’s \(t\) statistic is used as the “response”).

References

1. Adler Yu. P., Markova E. V., Granovsky Yu. V. Experimental Design in the Search for Optimal Conditions. Moscow, 1976.
2. Afifi A., Azen S. Statistical Analysis: A Computer-Oriented Approach. Translated from English. Moscow, 1982.
3. Bolch B., Huang K. J. Multivariate Statistical Methods for Economics. Translated from English. Moscow, 1979.
4. Voznesensky V. A. Statistical Methods of Experimental Design in Engineering and Economic Research. Moscow, 1981.
5. Gorbachev E. M., Kopanev V. A., Semenova V. N. In: Principles and Methods for Establishing Maximum Permissible Concentrations of Harmful Substances. Moscow, 1983, pp. 75–79.
6. Kovyazin V. G. Gig. truda, 1971, no. 2, pp. 54–56.
7. Kovyazin V. G. Gig. i san., 1978, no. 12, pp. 78–80.
8. Lyublina E. I., Minkina N. A., Rylova M. L. Adaptation to Industrial Poisons as a Phase of Intoxication. Leningrad, 1971.
9. Minchenko V. A. Gig. i san., 1984, no. 8, pp. 52–55.
10. Frolova A. D., Lisman M. B., Dvorkin E. A., Berliner E. G. Ibid., no. 10, pp. 56–58.
11. Shtabsky B. M., Kazan Yu. S., Katsnelson B. A. Ibid., 1983, no. 11, pp. 74–76.

Received June 3, 1986.