It is known that vascular endothelial growth factor (VEGF) expression is a response to hypoxia. On the other hand hypoxia may be detected by oximetry parameters including venous CO-oximetry indices or corresponding partial pressures of O_{2} and CO2. However significant correlation ties between VEGF levels and oximetry parameters were not found in groups of patients with ischemic stroke and transient ischemic attack. At that some effect related to the relationship between VEGF and sO_{2} was observed at corresponding scatter plots. Correlation between VEGF and proteins S100 levels in serum existed only in group with severe hypoxia where sO_{2} is less threshold close to 39-40%. So the relationship between VEGF level and saturation index sO_{2} exists in conjunction with additional factor that is S100 level in serum. To assess statistical significance of observed regularity it is necessary to test three null hypotheses about independence of one of involved factors on two another. The relationship may be manifestation of hypoxia effect on VEGF. To assess significance of hypoxia effect as a whole all three null hypotheses were tested with the help of developed technique based on random permutations and involving nonparametric combinations of criteria related to single oximetry parameters. The statistical significance assessments also involved multiplicity adjustment aimed to take into account multiple search of additional factors among variety of biological indices from analyzed data set. As a result of developed technique application all three considered null hypotheses were rejected at adjusted level p < 0.02 when effect of hypoxia on correlation between VEGF and complement component C4 was evaluated.
VEGF, Hypoxia, Controlled correlations, Permutation test, Nonparametric combinations, Multiple testing
Angiogenesis may be a biological response to insufficient oxygen supply resulting in hypoxia. The key mediator of angiogenesis and probably neurogenesis [1] is vascular endothelial growth factor (VEGF) which is homodimeric glycoprotein with a molecular weight of approximately 45 kDa. VEGF expression is activated as a response to stabilization and nuclear translocation of hypoxia-inducible factor-1 (HIF-1) when intracellular oxygen level is reduced [2]. Existence of HIF-1, VEGF signalling pathways is confirmed by high levels of VEGF in patients with chronic obstructive pulmonary disease (CHOPD) [3] and asthma [4] or in subjects with "plateau red face" [5]. It is known that immune system is involved in angiogenesis via secreting VEGF and other pro-angiogenic factors by macrophages, neutrophyles and other immune cells [6-8]. Oxygen saturation index sO_{2} and other CO-oximetry parameters in venous blood reflect balance between oxygen delivery and oxygen consumption [9]. Low sO_{2} correspond to tissue hypoxia [10]. However correlation coefficients between CO-oximetry parameters and VEGF levels were small and not statistically significant in groups of patients with severe neurological disorders. Corresponding data set is discussed below. Correlation coefficients between VEGF and partial pressures of O_{2} and CO_{2} were also not significant. Lack of reliable ties may be attributed among other things to complexity of existing dependence when relationship between two factors is controlled by the third one. Previous study of the data set with the help of OVP method [11,12] and visual analysis of related sparse diagrams uncover complex effect involving sO_{2} and levels in serum of VEGF and proteins from S100 family.
Left scatter diagram from Figure 1 corresponds to group with sO_{2} lower than 38.4%. This diagram conforms to existence of linear dependence of VEGF on S100 level. The only exception is one object that is marked with a red circle. At that there is no noticeable correlation between VEGF and S100 in the right diagram. This diagram corresponds to group with sO_{2} greater than 38.4%. It may be suggested that the correlation observed in the left diagram may be caused by severe hypoxia existing when sO_{2} is below 38.4%. The relationship from Figure 1 may be described by following equations:
Figure 1: Correlation between S100 and VEGF levels are compared in groups with sO_{2} < 38.4% (left scatter plot) and sO_{2} > 38.4% (right scatter plot). View Figure 1
$$\begin{array}{l}{y}_{j}\text{}=\text{}{\beta}_{0}^{l}\text{}+\text{}{\beta}_{1}^{l}{x}_{ij}\text{}+\text{}{\in}_{j}^{l}\text{}i\text{}f\text{}{z}_{kj}\text{}\text{}\delta \text{(1)}\\ {y}_{j}\text{}=\text{}{\beta}_{0}^{r}\text{}+\text{}{\beta}_{1}^{r}{x}_{ij}\text{}+\text{}{\in}_{j}^{r}\text{}i\text{}f\text{}{z}_{kj}\text{}\ge \text{}\delta \end{array}$$
Model (1) evidently is equivalent to standard piece-wise regression if X_{i} is equal to Z_{k}. It is supposed that the verification procedure must satisfy following demands:
• Verification procedure must include testing significance of both variables X and Z.
• Hypoxia is assessed by several oximetry parameters. It may be expected that statistical technique evaluating the effect of hypoxia on relationship between VEGF and additional factor X_{i} would be more powerful if it implements combination of statistical tests assessing effect each of oximetry parameters on relationship between Y and X_{i}.
• Effect from Figure 1 was found via testing variety of factors. So multiple testing must be taken into account when statistical significance is estimated.
Permutation test is an approach capable of meeting the listed demands due to following advantages. Permutation test may be implemented regardless whether or not underlying distributions of test statistics are known. There are no limitations on data sets sizes. Different permutation tests were proposed for assessing importance of each predictor in multiple regression models. Method based on testing significance of corresponding partial correlation coefficients may be mentioned in this regard [13]. Permutation test was used to make a choice between piece-wise regression and a simple linear regression [14]. However mentioned techniques are not suitable to evaluated significance of the effect considered in the paper. Verification is possible with the help of discussed below procedure based on testing several null hypotheses.
Effective methods calculating nonparametric combination (NPC) of several dependent permutation tests [15] may be used for assessing significance of hypoxia effect by several oximetry parameters. Also permutation tests are widely used to control multiplicity in various applied tasks including high-dimensional tasks related to DNA microarray experiments [16-19]. At that there are different ways to control FWER. Single step or stepwise procedures may be used to receive adjusted p-values from previously received raw p-values by min P correcting procedure or to calculate adjusted p-values directly from distributions of test statistics by maxT technique [20-22]. Unlike the mentioned works our paper is focused at problems associated with multiplicity control and combinations of several criteria when more complicated multifactor regularities are studied.
Effect of hypoxia on relationship between VEGF and 138 different biochemical, clinical or biophysical parameters was studied in a group of 88 patients of age from 33 to 88 with ischemic stroke and transient ischemic attack. Hypoxia level was assessed by partial pressures of O_{2} and CO_{2} and also by CO-oximetry parameters in venous blood that were measured with ABL80 FLEX CO-OX analyser. Serum levels of VEGF, S100 and complement component C4 were measured by the enzyme-linked immunosorbent assay (ELISA). CO-oximetry parameters together with partial pressures of O_{2} and CO_{2} in venous blood will be for the simplicity further referred to as oximetry parameters.
Results of standard correlation analysis are given in Table 1.
Table 1: Correlation coefficients between partial pressures of O_{2} and CO_{2} and CO-oximetry parameters in venous blood and VEGF levels. View Table 1
It can be seen from Table 1 that statistically significant linear ties between VEGF levels and hypoxia are absent.
It is quite possible that linear ties nevertheless may exist in combination with some additional factors. Such supposition is supported by the effect of oxygen saturation index (sO2) on relationship between VEGF and proteins S100 levels in serum. This effect was discovered with the help of OVP technique [11] and visual analysis of corresponding scatter plots. Left scatter plot from Figure 1 presents strong linear relationship between VEGF and S100 in group with saturation level sO_{2} less than 38.4. It is seen from right scatter plot that noticeable linear dependence between VEGF and S100 is absent in the group with saturation level sO_{2} greater than 38.4.
It may be supposed from Figure 1 that hypoxia leads to correlation between VEGF (Y ) and S100 (X_{i}). So S100 may be considered as additional factor. Our goal is to assess statistical significance of the assumed effect.
Complex dependencies testing: The above supposition is too complex to test using a single null hypothesis. In fact the effect contradicts the following three hypotheses:
a) ${H}_{0}^{1}(k,i)$ -Yis independent of vector (Z_{k}, X_{i});
b) ${H}_{0}^{2}(k,i)$ - X_{i} is independent of vector (Y, Z_{k});
c) ${H}_{0}^{3}(k,i)$ - Z_{k} is independent of vector (Y, X_{i}).
All three hypotheses must be rejected to be sure that the supposition is perfectly correct. Otherwise the observed pattern may be explained simpler. For example it may be attributed to existence of linear relationship between Y and X_{i} when ${H}_{0}^{3}(k,i)$ is in fact true.
Testing single hypothesis: The discussed effect is associated with great difference between ${\rho}_{l}\left(Y,{X}_{i}\right)$ and ${\rho}_{r}\left(Y,{X}_{i}\right)$ great value of one of the two correlation coefficients. High correlation values are achieved randomly in small groups. So the two conditions above may be a convincing evidence of the effect only when sizes of both groups are great enough. Thus the mentioned conditions correspond to great values of the following functional:
$${Q}_{2\rho}\left(\delta \right)\text{=}\frac{{m}_{l}{m}_{r}\left|\left|{\rho}_{l}\left(Y,{X}_{i}\right)\right|\text{}-\text{}\left|{\rho}_{r}\left(Y,{X}_{i}\right)\right|\text{}\right|}{1\text{}-\text{max}\left\{{\rho}_{l}^{2}\left(Y,{X}_{i}\right),\text{\hspace{0.33em}}{\rho}_{r}^{2}\left(Y,{X}_{i}\right)\right\}},\text{(2)}$$
Where m_{l} is the number of patients in group $\tilde{s}$_{l} with Z_{k} < δ and mr is the number of patients in group $\tilde{s}$r with Z_{k} > δ. The threshold δ is initially unknown. It is proposed to use the optimal threshold δ_{0} that is calculated as δ_{0} = arg max $\widehat{Q}$_{2ρ}(δ). Maximum of the functional Q_{2ρ}(δ) may be searched by trying all boundaries between distinct values of Z_{k} existing in full group $\tilde{s}$. Great values of Q_{2ρ}(δ_{0}) better testify against each of 3 null hypotheses ${H}_{0i}^{1}(k,i)$, ${H}_{0}^{2}(k,i)$ and ${H}_{0}^{3}(k,i)$ if δ_{0} is searched inside interval (δ_{1}; δr) including only such thresholds for Z_{k} that provide simultaneous validity of inequalities ml > mthr and mr > mthr. The narrowed search interval is used because of the high probability of great Q_{2ρ}(δ_{0}) values when null hypotheses are true but one of the groups is small. Probability of great correlation coefficients in small groups evidently is higher. Full compensation of this probability increase is impossible by using multiplier m_{l}m_{r} only.
Existence of outlying observations such as observation circled in red on the left part of Figure 1 reduces Q_{2ρ} and so may hinder correct statistical evaluation of the studied effect. So it is better to use robust Pearson correlation coefficient.
Robust correlation coefficient: Suppose that dependence of variable U on variable V is studied on data set $\tilde{s}$ = {(u_{j}, v_{j}) | j = 1, . . . , m}. At the first step simple linear regression model U = β_{0} + β_{1}V + e is calculated by $\tilde{s}$. Let $\tilde{s}$_{out} be the set of all 3σ outliers or $\tilde{s}$_{out} = {(u_{j}, v_{j}) || u_{j} - β_{0} - β_{1}v_{j}| > 3σ_{e} }, where σ_{e} is standard error of regression model. Robust correlation coefficient $\widehat{\rho}$ is calculated as standard Pearson correlation coefficient ρ by set $\tilde{s}$ \ $\tilde{s}$_{out}.
The robust version of the functional 2 can be written as
$${\widehat{Q}}_{2\rho}\left(\delta \right)\text{}=\text{}\frac{{m}_{l}{m}_{r}\left|\left|{\widehat{\rho}}_{l}\left(Y,{X}_{i}\right)\right|\text{}-\text{}\left|{\widehat{\rho}}_{r}\left(Y,{X}_{i}\right)\right|\right|}{1\text{}-\text{}\mathrm{max}\text{}\left\{{\widehat{\rho}}_{l}^{2}\left(Y,{X}_{i}\right),\text{}{\widehat{\rho}}_{r}^{2}\left(Y,{X}_{i}\right)\right\}},\text{(3)}$$
Where ${\widehat{\rho}}_{l}$ and ${\widehat{\rho}}_{r}$ are robust correlation coefficients in groups ${\tilde{s}}_{l}$ and $${\tilde{s}}_{r}$$. Optimal threshold δ_{0} is calculated now as ${\delta}_{0}\text{}=\text{}\mathrm{arg}{\mathrm{max}}_{\delta \in \left({\delta}_{l},{\delta}_{r}\right)}{\widehat{Q}}_{2\rho}\left(\delta \right)$. Functional ${\widehat{Q}}_{2\rho}\left({\delta}_{0}\right)$ is used as statistics for testing null hypotheses ${H}_{0}^{1}\left(k,i\right)$, ${H}_{0}^{2}\left(k,i\right)$ and ${H}_{0}^{3}\left(k,i\right).$ Then p-values for these null hypotheses are calculated according to the following Procedure I:
• Calculate optimal threshold δ_{0} on data set ${\tilde{S}}_{t}$ as $\mathrm{arg}{\mathrm{max}}_{\left({\delta}_{l},{\delta}_{r}\right)}{\widehat{Q}}_{2\rho}\left(\delta \right).$ The observed statistics value T0 is taken equal to ${\widehat{Q}}_{2\rho}\left({\delta}_{0}\right).$
To test ${H}_{0}^{1}\left(k,i\right)$ repeat independently steps a_{1} and a_{2} for r = 1, . . . , N.
• a_{1}) Take a random permutation ${\tilde{S}}_{t}^{r}$ of ${\tilde{S}}_{t}$. This is obtained by considering a random permutation u^{r} = $\left({u}_{1}^{r},\mathrm{...},{u}_{m}^{r}\right)$ of objects labels (1, . . . , m).
Let ${\tilde{S}}_{t}^{r}\text{=}\left\{\left({y}_{{u}_{1}^{r}},\text{}{x}_{1i},\text{}{z}_{1k}\right),\mathrm{...},\left({y}_{{u}_{m}^{r}},\text{}{x}_{mi},\text{}{z}_{mk}\right)\right\}.$
• a_{2}) Calculate optimal threshold δ_{0} on training set ${\tilde{S}}_{t}^{r}$ as $\mathrm{arg}{\mathrm{max}}_{\delta \in \left({\delta}_{l},{\delta}_{r}\right)}{\widehat{Q}}_{2\rho}\left(\delta \right).$ The statistics value T^{r} is taken equal to ${\widehat{Q}}_{2\rho}\left({\delta}_{0}\right).$
• Then calculate estimate of p-value for the null ${H}_{0}^{1}\left(k,i\right)$ as ${\widehat{p}}^{1}\left(k,i\right)\text{=}\frac{{\sum}_{1\le r\le N}\parallel \left({T}^{r}\ge {T}^{0}\right)}{N}.$
To ${H}_{0}^{2}\left(k,i\right)$ repeat independently steps a_{3} and a_{2} for r = 1, . . . , N.
• a3) Take a random permutation ${\tilde{S}}_{t}^{r}$ of ${\tilde{S}}_{t}^{}$ as ${\tilde{S}}_{t}^{r}\text{=}\left\{\left({y}_{1},{x}_{{u}_{1}^{r}i},{z}_{1k}\right),\mathrm{...},\left({y}_{m},{x}_{{u}_{m}^{r}i},{z}_{mk}\right)\right\}$ Where u^{r} = $\left({u}_{1}^{r},\mathrm{...},{u}_{m}^{r}\right)$ is a random permutation of {1, . . . , m}.
• Implement step (a_{2}).
• Then calculate estimate of p-value for the null ${H}_{0}^{2}\left(k,i\right)$ as ${\widehat{p}}^{2}\left(k,i\right)\text{=}\frac{{\sum}_{1rN}\parallel \left({T}^{r}{T}^{0}\right)}{N}.$
To test ${H}_{0}^{3}\left(k,i\right)$ repeat independently steps a_{4} and a_{2} for r = 1, . . . , N.
• a_{4}) Take a random permutation ${\tilde{S}}_{t}^{r}$ of ${\tilde{S}}_{t}^{}$ as ${\tilde{S}}_{t}^{r}\text{=}\left\{\left({y}_{1},{x}_{11},{z}_{{u}_{1}^{r}}\right),\mathrm{...},\left({y}_{m},{x}_{m1},{z}_{{u}_{m}^{r}}\right)\right\}$ Where u^{r} = $\left({u}_{1}^{r},\mathrm{...},{u}_{m}^{r}\right)$ is a random permutation of {1, . . . , m}.
• Implement step (a_{2}).
• Then calculate estimate of p-value for null hypothesis ${H}_{0}^{3}\left(k,i\right)$ as ${\widehat{p}}^{3}\left(k,i\right)\text{=}\frac{{\sum}_{1rN}\parallel \left({T}^{r}{T}^{0}\right)}{N}.$
A regularity for combination of factors (Y, X_{k}, Z_{i}) similar to the one in Figure 1 is supposed to be significant at level α if all three inequalities ${\widehat{p}}^{1}\left(k,i\right)<\alpha ,{\widehat{p}}^{2}\left(k,i\right)<\alpha \text{and}{\widehat{p}}^{3}\left(k,i\right)\alpha $ hold. It is difficult to make a theoretical conclusion about unbiasedness and consistency of the described test. However its performance can be assessed in experiments with simulated data.
Design of experiments: Experiments were designed to imitate regularity from Figure 1 only for certain groups of variables while for the remaining ones such regularities were absent. Scenario includes generating variables Y, Z_{1}, . . . , Z_{7}, X_{1}, . . . , X_{n}. VariablesYand e_{1}, . . . , e_{n} were independently sampled from N(0,1), variable X_{1}, U_{g} and variables U_{2}, ..., U_{7} are independently sampled from continuous uniform distribution U(0,1). Variables X_{1}, . . . , X_{n} are calculated fromY, e_{1}, . . . , e_{n}, Z_{1} and U_{g}. Variables Z_{2}, . . . , Z_{7} are calculated from Z_{1} and U_{2}, . . . , U_{7}.
1) Variables X_{1}, . . . , X_{30} for an observation j were calculated as ${X}_{ij}\text{=}{Y}_{j}-{\delta}_{1}\text{}\ast \text{}{e}_{ij}$ if Z_{1j} < 0.4 and as ${X}_{ij}\text{=}{\delta}_{2}\text{}\ast \text{}{e}_{ij}$ otherwise. In groups with Z_{1j} < 0.4 correlation level betweenYand X_{i} is determined by parameter δ_{1}. Experiments with δ_{1} = 0.75 and δ_{1} = 0.57 were conducted. Choice δ_{1} = 0.75 provides generating data with ρ(Y, X_{i}) = 0.8 and choice δ_{1} = 0.57 corresponds to ρ(Y, X_{i}) = 0.87. Parameter δ_{2} was taken equal to $\sqrt{1\text{+}{\delta}_{1}^{2}}.$ Observations with Z_{1j} ≥ 0.4 were generated from distribution with ρ(Y, X_{i}) = 0.
Thus in the first experiment data was sampled from distribution with ρ(Y, X_{i}) = 0.8, σ_{Y} = 1 and ${\sigma}_{{X}_{i}}\text{=}\sqrt{1\text{+}{0.75}^{2}}$ when Z_{1} < 0.4 and from distribution with ρ(Y, X_{i}) = 0, σ_{Y} = 1 and ${\sigma}_{{X}_{i}}\text{=}\sqrt{1\text{+}{0.75}^{2}}$ when Z_{1} ≥ 0.4.
In the second experiment data was sampled from distribution with ρ(Y, X_{i}) = 0.87, σ_{Y} = 1 and ${\sigma}_{{X}_{i}}\text{=}\sqrt{1\text{+}{0.57}^{2}}$ when Z_{1} < 0.4 and from distribution with ρ(Y, X_{i}) = 0, σ_{Y} = 1 and ${\sigma}_{{X}_{i}}\text{=}\sqrt{1\text{+}{0.57}^{2}}$ when Z_{1} ≥ 0.4. So data is generated to provide existence of effect that is similar to effect from Figure 1 for each combination from set {(Y, X_{i},Z_{k})|i = 1, . . . , 30, k = 1, . . . , 7}.
2) Variables X_{31}, . . . , X_{60} for the observation j were calculated as ${X}_{ij}\text{=}{Y}_{j}-{\delta}_{1}\ast {e}_{ij}$ if U_{g} < 0.4 and as ${X}_{ij}\text{=}{\delta}_{2}\ast {e}_{ij}$ if U_{g} ≥ 0.4. So observations are generated from mixture of distributions with ρ(Y, X_{i}) = 0.8 and ρ(Y, X_{i}) = 0 in the first experiment and from mixture of distributions with ρ(Y, X_{i}) = 0:87 and ρ(Y, X_{i}) = 0 in the second one. Thus data is generated to provide existence of weak mutual correlation for pairs from {(Y, X_{i})|i = 31, . . . , 60} and independence of such pairs on Z_{1}, . . . , Z_{k}.
3) Variables X_{61}, . . . , X_{90} for the observation j were calculated as ${X}_{ij}\text{=}{\delta}_{2}\ast {e}_{ij}.$ SoYis independent on variables X_{61}, . . . , X_{90} and Z_{1}.
4) Variables Z_{2}, . . . , Z_{7} for the observation j were calculated as Z_{ij} = Z_{1j} if U_{g} ≤ 0.9 and as Z_{ij} = 1 - Z_{1j} if U_{g} > 0.9. There are no regularities for combinations (Y,X_{k},Z_{i}) when 31 ≤ i ≤ 90 that are similar to regularity from Figure 1. Variables Z_{1}, . . . , Z_{7} are included in scenario to imitate all 7 oximetry factors. At that for combinations (Y,X_{1},Z_{i}) regularities are more distinct to compare with regularities for combinations (Y,X_{k},Z_{i}) when k > 1. Sets X_{1}, . . . , X_{30}; X_{31}, . . . , X_{60}; X_{61}, . . . , X_{90} will be referred to as ${\tilde{C}}_{1},\text{}{\tilde{C}}_{2},\text{}{\tilde{C}}_{3}$ correspondingly.
Results of experiment: Results of the first and second experiments are presented in Table 2. Columns of the table correspond to significance levels from p < 0.0001 to p < 0.1. Upper part of table corresponds to first experiment (δ_{1} = 0.75) and lower part of table corresponds to second experiment (δ_{1} = 0.57). Cell at intersection of row corresponding to significance level α and column corresponding to subset ${\tilde{C}}_{j}$ contains number of triples (Y, X_{i},Z_{1}) with ${X}_{i}\text{}\in \text{}{\tilde{C}}_{j}$ for which all three null hypotheses were rejected at least at level α. Number of triples from set $\left\{\left(Y,{X}_{i},{Z}_{k}\right)|{X}_{i}\text{}\in \text{}{\tilde{C}}_{j},2\le k\le 7\right\}$ for which all three null hypotheses were rejected at level α is given in the same cell in parentheses.
Table 2: Results of experiments with simulated data. View Table 2
It may be seen from Table 2 that all three null hypotheses were rejected for the majority combinations from the set ${\tilde{C}}_{1}\text{=}\left\{\left(Y,{X}_{i},{Z}_{k}\right)|1\text{}\le \text{}i\text{}\le \text{30,}k\text{=1,}\text{.}\text{.}\text{.,}7\right\}.$ In the first experiment all null hypotheses were rejected at significance level p < 0.01 for all 30 combinations from ${\tilde{C}}_{1}$ when k = 1 and for 163 combinations from 180 when k = 2, . . . ,7. On the contrary all null hypotheses were rejected only for few combinations from the sets ${\tilde{C}}_{2}\text{=}\left\{\left(Y,{X}_{i},{Z}_{k}\right)|31\text{}\le \text{}i\text{}\le \text{60,}k\text{=1,}\text{.}\text{.}\text{.,}7\right\}$ and ${\tilde{C}}_{3}\text{=}\left\{\left(Y,{X}_{i},{Z}_{k}\right)|61\text{}\le \text{}i\text{}\le \text{90,}k\text{=1,}\text{.}\text{.}\text{.,}7\right\}.$ All null hypotheses were rejected at significance level p < 0.1 only for one combination from ${\tilde{C}}_{2}$ and for one combination from ${\tilde{C}}_{3}$ when k = 1. At that number of combinations where all three null hypotheses were rejected equals 6 inside ${\tilde{C}}_{2}$ and 3 inside ${\tilde{C}}_{3}$ for k = 2, . . . ,7. In the second experiment number of combinations where three null hypotheses were rejected at level α were higher than number of such combinations practically for all significance levels. So the results of experiments strongly indicate unbiasedness of the developed criterion.
The developed technique was applied to find regularities similar to the one shown in Figure 1 on the described above clinical data set. Three null hypotheses were tested for combinations (Y, X_{i},Z_{k}) whereYis concentration of VEGF in serum, Z_{1}, . . . , Z_{7} were oximetry parameters sO2, pO2, pCO2, FCOHb, FO2Hb, FMetHb, FHHb. All 138 variables different from VEGF concentration and oxymetry parameters were tried as additional factors X_{1}, . . . , X_{n}. The most significant effects were revealed if concentrations of S100 proteins or complement component C4 in serum are used as additional factor X. Table 3 present calculated p-values ${\widehat{p}}^{1},\text{}{\widehat{p}}^{2}\text{and}{\widehat{p}}^{3}$ correspondingly for all such combinations.
Table 3: Results of the null hypotheses testing when S100 or C4 are additional factors. View Table 3
It is seen from Table 3 that all null hypotheses are rejected at significance level p < 0.001 for combination (VEGF, pO2, C4), at significance level p < 0.002 for combinations (VEGF, sO2, C4) and (VEGF, FO2Hb, C4), at significance level p < 0.05 for combinations (VEGF, FHHb, C4), (VEGF, sO2, S100), (VEGF, pCO2, S100). Regularity related to effect of hypoxia on correlation between VEGF and S100 is shown at Figure 2.
Figure 2: Correlations between S100 and VEGF levels are compared in groups with sO_{2} < 39.75% (left scatter plot) and sO_{2} > 39.75% (right scatter plot). Boundary point was calculated using Procedure I. View Figure 2
Pattern from Figure 2 is similar to the pattern from Figure 1. However boundary point for pattern from Figure 2 is received by procedure I. This boundary differs from boundary for the pattern from Figure 1 that was calculated by OVP method. Correlation coefficient between VEGF and S100 in group of 33 cases with sO_{2} < 39.75% equals 0.64. Correlation coefficient increases to 0.88 after removing an outlying object highlighted in the left scatter diagram by red circle. No relationship between VEGF and S100 exists in group of 55 cases with sO_{2} < 39.75% as it may be seen from right diagram. Corresponding correlation coefficient equals 0.03.
It may be seen from Figure 3 and Figure 4 that effect of hypoxia on relationship between VEGF and C4 is similar to the effect of hypoxia on relationship between VEGF and S100. Correlation coefficient between VEGF and C4 is equal 0.47 in group of 31 cases with sO_{2} < 39.25% which corresponds to left scatter diagram from Figure 3. Correlation coefficient increases to 0.76 after removing of outlying object which is highlighted at the left scatter diagram by red circle. No significant relationship between VEGF and C4 exists in group of 57 cases with sO_{2} > 39.25% as it may be seen from the right diagram. Corresponding correlation coefficient equals -0.11. Correlation coefficient increases to 0.05 after removing the highlighted at right diagram outlier.
Figure 3: Correlation between VEGF and C4 levels are compared in groups with sO_{2} < 39.25% (left scatter plot) and sO_{2} > 39.25% (right scatter plot). Boundary point was calculated using Procedure I. View Figure 3
Figure 4: Correlation between VEGF and C4 levels are compared in groups with FHHB < 56.3% (left scatter plot) and FHHB > 56.3% (right scatter plot). Boundary point was calculated using Procedure I. View Figure 4
Low saturation index sO_{2} corresponds to high FHHB values. Strong correlation between VEGF and S100 exists when FHHB is greater than a certain threshold. At that correlation coefficient is close to zero when FHHb is lower than the threshold as can be seen in Figure 4. Correlation coefficient between VEGF and C4 in group of 36 cases with FHHB > 56.3% equals 0.45. Correlation coefficient increases to 0.73 after removing of outlier highlighted at right diagram. No significant relationship between VEGF and C4 exists in group of 52 cases with FHHB < 56.3% as it may be seen from left diagram. Corresponding correlation coefficient is equal -0.1. Correlation coefficient increases to 0.07 after removing of highlighted at left diagram outlier.
Our goal is testing if hypoxia has effect on VEGF production via controlling relationship between VEGF and some additional factor. Hypoxia effect is manifested via effects associated with different oximetry parameters. Existence of supposed hypoxia effect contradicts simultaneously to several of null hypotheses associated with different oximetry parameters. Hypoxia effect may be assessed by testing global null hypotheses ${H}_{0c}^{1}\left(i\right)\text{=}{\cap}_{k=1}^{7}{H}_{0}^{1}\left(k,i\right),$ ${H}_{0c}^{2}\text{=}{\cap}_{k=1}^{7}{H}_{0}^{2}\left(k,i\right)$ and ${H}_{0c}^{3}\text{=}{\cap}_{k=1}^{7}{H}_{0}^{3}\left(k,i\right)$ [15]. Last global hypotheses may be tested with the help of nonparametric combinations (NPC) methodology [15]. To test global hypothesis ${H}_{c0}^{g}\left(i\right)\text{}\left(g\text{=}1\text{,}\text{.}\text{.}\text{.,3}\right)$ on data set ${\tilde{S}}_{t}\text{=}\left\{\left({y}_{1},\text{}{x}_{1i},\text{}{z}_{11},\text{}\text{.}\text{.}\text{.,}{z}_{17}\right),\text{}\text{.}\text{.}\text{.,}\left({y}_{m},\text{}{x}_{mi},\text{}{z}_{m1},\text{}\text{.}\text{.}\text{.,}{z}_{m7}\right)\right\}$ following NPC procedure was used.
Procedure II: Repeat steps (b1), . . . , (b4) for k = 1, . . . , 7.
• b_{1}) Calculate T0 according to step a0 from Procedure I on data set $\left\{\left({y}_{1},{x}_{1i},{z}_{1k}\right),\text{}\text{.}\text{.}\text{.,}\left({y}_{m},{x}_{mi},{z}_{mk}\right)\right\}$
• b_{2}) For r = 1, . . . , N repeat step if is tested, g = 1, . . . , 3.
- $\left.{b}_{2}^{1}\right)$ Take random permutation ${\tilde{S}}_{t}^{r}\text{=}\left\{\left({y}_{{u}_{1}^{r}},{x}_{1i},{z}_{1k}\right),\text{}\text{.}\text{.}\text{.,}\left({y}_{{u}_{m}^{r}},{x}_{mi},{z}_{mk}\right)\right\}$ of ${\tilde{S}}_{t}$ according to step (a_{1}) from Procedure I and implement step a_{2} from Procedure I to calculate statistics Tr by ${\tilde{S}}_{t}^{r}.$
- $\left.{b}_{2}^{2}\right)$ Take random permutation ${\tilde{S}}_{t}^{r}\text{=}\left\{\left({y}_{1},{x}_{{u}_{1}^{r}i},{z}_{1k}\right),\text{}\text{.}\text{.}\text{.,}\left({y}_{m},{x}_{{u}_{m}^{r}i},{z}_{mk}\right)\right\}$ of ${\tilde{S}}_{t}$ according to step (a_{3}) from Procedure I and implement step a_{2} from Procedure I to calculate statistics Tr by ${\tilde{S}}_{t}^{r}.$
- $\left.{b}_{2}^{3}\right)$ Take random permutation ${\tilde{S}}_{t}^{r}\text{=}\left\{\left({y}_{1},{x}_{1i},{z}_{{u}_{1}^{r}k}\right),\text{}\text{.}\text{.}\text{.,}\left({y}_{m},{x}_{mi},{z}_{{u}_{m}^{r}k}\right)\right\}$ of ${\tilde{S}}_{t}$ according to step (a_{4}) from Procedure I and implement step a_{2} from Procedure I to calculate statistics Tr by ${\tilde{S}}_{t}^{r}.$
• b_{3}) Calculate $\widehat{p}\left(k,i\right)\text{=}\frac{{\sum}_{1\le j\le N}\parallel \left({T}^{j}\left(i\right)\ge {T}^{0}\left(i\right)\right)}{N}.$
• b_{4}) Calculate $\lambda \left(k,r,i\right)\text{=}\frac{{\sum}_{1\le j\le N}\parallel \left[{T}^{j}\left(i\right)\ge {T}^{r}\left(i\right)\right]}{N}$ r = 1, . . . , N
• Calculate statistics ${T}_{c}^{0}\left(i\right)\text{=}\psi \left[\widehat{p}\left(1,i\right),\text{}\text{.}\text{.}\text{.,}\widehat{p}\left(7,i\right)\right]$
• Calculate statistics ${T}_{c}^{r}\left(i\right)\text{=}\psi \left[\lambda \left(1,r,i\right),\text{}\text{.}\text{.}\text{.,}\lambda \left(7,r,i\right)\right]$
• Then calculate ${\widehat{p}}_{c}\left(i\right)\text{=}\frac{{\sum}_{1\le j\le N}\parallel \left[{T}_{c}^{r}\left(i\right)\ge {T}_{c}^{0}\left(i\right)\right]}{N}$ that is p-value testing hypothesis ${H}_{0c}^{g}$ when step ${b}_{2}^{g}.$
Test statistics in Procedure II is calculated as a combining function of p-values related to partial tests. Several combining functions are discussed in [15]. The best performance is achieved according to our experiments when slightly modified Fisher combining function ψ is used. Let ${\widehat{p}}_{1},\text{}\text{.}\text{.}\text{.,}{\widehat{p}}_{k}$ are some p-values calculated by permutation test with N random permutations. Then $\psi \left({\widehat{p}}_{1},\text{}\text{.}\text{.}\text{.,}{\widehat{p}}_{k}\right)\text{=}-2{\sum}_{i=1}^{k}\mathrm{log}{\widehat{p}}^{\prime}{}_{i},$ where ${\widehat{p}}^{\prime}{}_{i}\text{=}{\widehat{p}}_{i}$ if ${\widehat{p}}_{i}\text{}0$ and ${\widehat{p}}^{\prime}{}_{i}\text{=1/2}N$ otherwise.
Results of the Procedure II applied to the studied data set are presented in Table 4. It can be seen that the global null hypotheses ${H}_{0c}^{1}$ and ${H}_{0c}^{2}$ are rejected at level p < 0.0005 when additional factors are C4 and S100 concentrations. Global null hypothesis ${H}_{0c}^{3}$ is rejected at p = 0.0001 when additional factor is C4 concentration. But ${H}_{0c}^{3}$ is not rejected when additional factor is S100 concentration.
It was necessary to test global null hypotheses from set $\left\{{H}_{0c}^{1},\text{}{H}_{0c}^{2},\text{}{H}_{0c}^{3}\left|i\text{=}1,\text{}\text{.}\text{.}\text{.,}n\right.\right\}$ to reveal the supposed effect of hypoxia on relationship between concentration of VEGF and concentrations of C4 or S100. So a multiple testing procedure must be used to assess true statistical significance of revealed effects. It is sufficient to use a single-step procedure because only two global null hypotheses associated with C4 and S100 as initial factors were rejected. Methods min P and max T are effective tools of multiplicity control [22]. More universal min P procedure is preferable for discussed task because max T technique is based on unjustified supposition about approximate equality of test statistics distribution.
Adjusted p-values were calculated by the represented below Procedure III.
Repeat steps (c_{1}), (c_{2}) and (c_{3}) for i = 1, . . . , n.
• c_{1}) Use Procedure II to calculate statistics ${T}_{c}^{0}\left(i\right),\text{}{T}_{c}^{1},\text{}\text{.}\text{.}\text{.,}{T}_{c}^{N}$ by ${\tilde{S}}_{t}$ and random permutations $\left\{{\tilde{S}}_{t}^{r}\left|r\text{}=\text{}1,\text{}.\text{}.\text{}.\text{},\text{}N\right.\right\}$ correspondingly. At that step ${b}_{2}^{g}$ is implemented when ${H}_{0c}^{g}\left(i\right)$ is tested, g = 1, . . . , 3.
• c_{2}) Calculate ${\lambda}_{c}\left(0,i\right)\text{=}\frac{{\sum}_{1\le j\le N}\parallel \left({T}_{c}^{j}\left(i\right)\ge {T}_{c}^{0}\left(i\right)\right)}{N}$
• c_{3}) Calculate ${\lambda}_{c}\left(r,i\right)\text{=}\frac{{\sum}_{1\le j\le N}\parallel \left({T}_{c}^{j}\left(i\right)\ge {T}_{c}^{r}\left(i\right)\right)}{N}$
• c_{4}) Calculate ${\widehat{p}}_{mt}\text{=}\frac{{\sum}_{1\le r\le N}\parallel \left[{\mathrm{min}}_{i\in \left\{1,\mathrm{...},n\right\}}{\lambda}_{c}^{r}\left(i\right)\le {\widehat{p}}_{c}^{0}\left(i\right)\right]}{N}$
• c_{5}) ${\widehat{p}}_{mt}$ is equal to ${\widehat{p}}_{mt}^{g}$ that is adjusted p-value testing global null hypothesis ${H}_{0c}^{g}$ when $\left({b}_{2}^{g}\right)$ is used at step (c_{1}).
The described Procedure III provides weak FWE control. Procedure III was applied to calculate adjusted p-values for global null hypotheses ${H}_{0c}^{1},$ ${H}_{0c}^{2}$ and ${H}_{0c}^{3}.$ Results are presented in Table 4.
Table 4: Results of the global null hypothesis testing with NPC and multiplicity control with Procedure III. View Table 4
It is seen from Table 4 that the global null hypotheses ${H}_{0c}^{1},$ and ${H}_{0c}^{2}$ are rejected when concentrations of C4 and S100 are additional factors. At that ${H}_{0c}^{3}$ is not rejected if additional factor is concentration of S100. So it is possible that pattern from the Figure 2 may be related only to existence of linear relationship between VEGF and S100 levels in serum which is not controlled by hypoxia. On the other hand set of relationships that are represented at Figure 3 and Figure 4 or are mentioned in Table 3 where concentrations of C4 is additional factor cannot be explained with the help of some simpler effect. These relationships cannot be reduced to linear correlation between VEGF and C4. Also they cannot be explained by effect of hypoxia on VEGF only or on C4 only. So supposition that relationship between VEGF and C4 levels in serum is controlled by hypoxia is in accordance with data.
Results may be shortly summarized as follows. A method was developed which is aimed to discover relationships of the following type in data: Significant linear correlation between two factorsYand X_{i} exists only if third factor Z_{k} belongs to interval from one side of some threshold δ. At that from another side of δ Pearson correlation coefficient is close to zero.
It was suggested to consider such three-factor relationship as statistically significant when rejecting three null hypotheses: ${H}_{0c}^{1},$ ${H}_{0c}^{2}$ and ${H}_{0c}^{3}.$ Nonparametric permutation tests with statistics that is optimal value of special quality functional were used to test these hypotheses.
Performance of the method was evaluated in tasks with simulated data. Good concordance between found regularities and patterns provided by the experiment scenario is seen from Table 2.
The method was applied to test supposition that hypoxia control relationship between serum VEGF concentration and some factor from the analyzed clinical database. It was supposed that hypoxia is manifested by oximetry parameters. Three null hypotheses were rejected for set of triples (Y, X_{i},Z_{k}) whereYis VEGF concentration, X_{i} is some additional factor and Z_{k} is some of oximetry parameters.
Significance of hypoxia effect on correlation between VEGF level and additional factor X_{i} may be assessed as combined significance of effects related to different oximetry parameters. Combined significance was evaluated with the help of NPC method testing intersection of null hypotheses related to set of triples $\left\{\left(Y,{X}_{i},{Z}_{k}\right)\left|k\text{=1,}\text{.}\text{.}\text{.,7}\right.\right\}.$ A single-step permutations based FWE control was implemented to take into account that additional factor is searched among 138 variables.
It was shown that three combined null hypotheses were rejected at significance level p < 0.02 when concentration of complement C4 is the additional factor. Developed technique may be used in variety of biomedical tasks where it is necessary to assess effect of some factor or some group of factors on existing linear ties.
This study was supported by RFBR grant 20-01-00609.