In genetic association studies, important and common goals are the
identification of single nucleotide polymorphisms (SNPs) showing a
distribution that differs between several groups and the detection of
SNPs with a coherent pattern. In the former situation, tens of thousands
of SNPs should be tested, whereas in the latter case typically
several ten SNPs are considered leading to thousands of statistics that
need to be computed.
A test statistic appropriate for both goals is Pearson’s Chi^2-statistic.
However, computing this (or another) statistic for each SNP or pair
of SNPs separately is very time-consuming.
In this article, we show how simple matrix computation can be
employed to calculate the Chi^2-statistic for all SNPs simultaneously

Schwender, Holger

Eldorado - Ressourcen aus und für Lehre, Studium und Forschung

A Note on the Simultaneous Computation of
Thousands of Pearson’s χ2−Statistics
Holger Schwender
Collaborative Research Center 475
Department of Statistics
University of Dortmund
holger.schwender@udo.edu
Abstract
In genetic association studies, important and common goals are the
identification of single nucleotide polymorphisms (SNPs) showing a
distribution that differs between several groups and the detection of
SNPs with a coherent pattern. In the former situation, tens of thou-
sands of SNPs should be tested, whereas in the latter case typically
several ten SNPs are considered leading to thousands of statistics that
need to be computed.
A test statistic appropriate for both goals is Pearson’s χ2-statistic.
However, computing this (or another) statistic for each SNP or pair
of SNPs separately is very time-consuming.
In this article, we show how simple matrix computation can be
employed to calculate the χ2-statistic for all SNPs simultaneously.
1
1 Introduction
While association studies typically comprise the genotypes of several ten sin-
gle nucleotide polymorphisms (SNPs), quite recently developed microarrays
allow to measure the genotypes of tens or even hundreds of thousands of
SNPs simultaneously. In the former situation, two tasks are the identifica-
tion of SNPs showing a distribution that differs substantially between several
groups (e.g., non-cancer vs. cancer) and the detection of groups of SNPs with
a coherent pattern. Considering SNP microarray, a first goal is to reduce the
number of SNPs to a better manageable size.
Since the latter goal is similar to the first task in the former situation,
both problems can be solved in the same way: For each SNP, a statistic
appropriate for testing if its distribution differs between several groups is
computed. The higher this score, the more likely it is that the corresponding
SNP differs substantially between the classes (for topics such as adjusting for
multiplicity that have to be considered in this situation, see, e.g., Dudoit et
al., 2003).
As SNPs are categorical variables exhibiting three realizations (homozy-
gous reference, heterozygous, and homozygous variant), an appropriate test
statistic is Pearson’s χ2-statistic
χ2 =
r∑
g=1
c∑
k=1
(ngk − n˜gk)2
n˜gk
, (1.1)
where ngk and n˜gk are the observed and expected numbers of observations,
respectively, shown in the gth row and kth column of the corresponding con-
tingency table, g = 1, . . . , r, k = 1, . . . , c. Here, c = 3, and r is the number
of groups such that ngk specifies how many of the n =
∑
g,k ngk =
∑
g,k n˜gk
observations in the gth class showing the kth genotype.
A solution to the other problem, i.e. the detection of groups composed of
2
related SNPs, is, again, to compute Pearson’s χ2-statistic (1.1) with r = c = 3
to test if two SNPs are independent. Afterwards, e.g., Pearson’s corrected
contingency coefficient
PC =
√
d
d− 1 ·
χ2
χ2 + n
with d = min{r, c} here= 3 is calculated as a measure of similarity for these
two SNPs.
Since considering each pair of m SNPs means that m(m− 1)/2 similari-
ties/distances have to be determined, Pearson’s χ2-statistic has to be calcu-
lated several hundred to a few thousand times even if m < 100. The same
applies to the analysis of microarrays in which tens of thousands of SNPs
are tested for different group distributions. In these cases, it can therefore
be time-consuming to compute each of the χ2-statistics separately.
In this article, we show how simple matrix calculation can be employed
to consider all SNPs or pairs of SNPs simultaneously. Another advantage of
this approach is that it provides a matrix composed of all m or m(m− 1)/2
contingency tables, respectively. This matrix thus enables the determination
of other similarity measures based on contingency tables such as simple or
flexible matching coefficients (Mu¨ller et al., 2005).
This article is organized as follows: In Section 2, the approaches for the
simultaneous computation of thousands of Pearson’s χ2-statistics in the two
situations are described, whereas in Section 3 we discuss practical issues
as the actual implementation of these algorithms and the handling of both
missing values and variables with differing numbers of levels. This section
also contains a short description of how the approach for detecting groups of
SNPs with coherent patterns can be used to compute similarity measures that
3
are not based on Pearson’s χ2-statistic. Finally, the matrix based algorithms
are compared with the corresponding individual computations in Section 4.
2 Simultaneous Computation of Pearson’s χ2-
Statistic
Let X be an m × n matrix in which each column corresponds to one of the
n observations and each row to one of the m variables, and y be a vector of
length n containing the class labels of the n observations. Assume that each
of these variables exhibits c levels denoted by the integers 1, . . . , c, and that
each observation belongs to one of the r classes 1, . . . , r. Then, Pearson’s
χ2-statistics for testing a variable if its distribution differs between r groups
can be computed for all variables represented in X simultaneously by the
procedure described in Algorithm 1.
Algorithm 1 (Rowwise Pearson’s χ2-Statistic)
Let X be an m × n matrix consisting of the values 1, . . . , c, and y be a
vector of length n containing the class labels 1, . . . , r of the n observations
represented by the columns of X.
1. Let X(k) denote an m× n matrix with elements
x
(k)
ij =
1, if xij = k0 otherwise ,
k = 1, . . . , c, and L be an n× r matrix with elements `jg = I
(
yj = g
)
.
2. For k = 1, . . . c, set L(k) = X(k)L and
L˜(k) =
1
n
X(k)1nL
′1n, (2.1)
4
where 1n is a vector consisting of n ones, and compute
S(k) =
L(k) ∗ L(k)
L˜(k)
(2.2)
with s
(k)
ig = `
(k)
ig · `(k)ig / l˜(k)ig , i = 1, . . . ,m, g = 1, . . . , r.
3. Let S be an n× c matrix in which the kth column consists of the vector
sk = S
(k)1r, k = 1, . . . , c. The vector comprising Pearson’s χ
2-statistics
for testing each row of X if the distribution of the corresponding vari-
able differs between the groups specified by y is given by
rX = S1c − n. (2.3)
For the computations in Algorithm 1, note that (1.1) can also be expressed
as
χ2 =
r∑
g=1
c∑
k=1
n2gk
n˜gk
− n,
and that the
(
ith, gth
)
element of L(k) or L˜(k) comprises the observed or
expected number of observations, respectively, being a member of group g
and showing the kth level at the ith variable.
Using similar ideas, Algorithm 2 describes how Pearson’s χ2-statistic for
testing if two variables are independent can be determined for all of the
m(m − 1)/2 pairs of m variables simultaneously. Note that Algorithm 2
assumes that all variables exhibit the same number of levels such that r = c.
For an extension to r ≤ c, see Section 3.3.
5
Algorithm 2 (Pairwise Pearson’s χ2-Test for Independence)
Let X be an m× n matrix consisting of the values 1, . . . , c.
1. Let X(k) denote an m× n matrix, k = 1, . . . , c, with elements
x
(k)
ij =
1, if xij = k0 otherwise .
2. For g, k = 1, . . . , c, compute N(gk) = X(g)X(k)′ and
N˜(gk) =
1
n
X(g)1n1
′
nX
(k)′. (2.4)
3. Pearson’s χ2-statistic for testing if the ith and the hth variable, i, h =
1, . . . ,m, represented in the ith and hth row of X, respectively, are
independent is given by the
(
ith, hth
)
element of
RX =
c∑
g=1
c∑
k=1
N(gk) ∗N(gk)
N˜(gk)
− n. (2.5)
3 Practical Issues
3.1 Details on the Actual Implementation
Algorithms 1 and 2 have been implemented in the statistical software envi-
ronment R (Ihaka and Gentleman, 1996). The following notes give details on
their actual implementations:
6
• Instead of computing the rowwise or columnwise sums of a matrix by
multiplying it with a vector consisting of an appropriate number of
ones, the even faster R function rowSums or colSums are employed for
these calculations.
• In (2.2) and (2.5), the elementwise squaring of L(k) or N(gk) is repre-
sented by “*”, since in R it is faster to elementwise square a matrix Z
by Z ∗ Z than by Z2. Thus, the actual reason for this is not to avoid
notations such as N(gk)2.
• In the second step of Algorithm 2, not all c2 matrices N(gk) are deter-
mined. We only consider N(gk) for g = 1, . . . , c, and k = g, . . . , c, as the
lower (upper) triangle of N(gk) contains the same values as the upper
(lower) triangle of N(kg). The same applies to N˜(gk).
• Finally, the pairwise χ2-statistics are actually not computed as shown
in (2.5). Instead, the upper triangle of N(gk) is stored in the
(
(g −
1)c + k
)th
column of the m(m − 1)/2 × c2 matrix M, and the lower
triangle in the
(
(k− 1)c+ g)th column. Hence, this matrix M contains
all contingency tables corresponding to any of the m(m − 1)/2 pairs
of m variables. Therefore, M can also be employed to compute other
similarity measures based on contingency tables. In the same way, all
expected cell entries are retained in a matrix M˜.
The vector rX comprising the m(m − 1)/2 χ2-statistics can thus be
determined by
rX =
(
M ∗M
M˜
− n
)
1c2 .
By default, rX is then stored in the lower triangle of an m×m matrix
such that the lower triangle of this matrix is identical to the lower
triangle of RX in (2.5).
7
3.2 Missing Values
Even though both Algorithm 1 and 2 still work if there are missing values,
the resulting χ2-statistics will not be correct, since the expected numbers of
observations are divided by n, i.e. the total number of observations, and not
by the actual number of observations showing no missing value at a particular
variable.
For solving this problem, let XA be an m× n matrix with
xAij =
1, xij is not missing0, xij is a missing value ,
and replace n in (2.3) by the vector
n =

n1
...
nm
 = XA1n.
In (2.1), 1
n
is substituted by the m× r matrix
Nden =

1
n1
· · · 1
n1
...
. . .
...
1
nm
· · · 1
nm
 ,
whereas in (2.4) and in (2.5) n is replaced by the m×m matrix
N = XA
(
XA
)′
.
Since (2.4) just takes individual and not pairwise missing values, i.e. miss-
ing values appearing in either of the two considered variables, into account,
it is additionally necessary to replace the rowwise sums of X(k) by Z(k) =
X(k)
(
XA
)′
such that (2.4) becomes
N˜(gk) =
Z(g) ∗ (Z(k))′
N
.
8
3.3 Different Numbers of Levels
If the ith variable exhibits r < c levels, where c is the maximum number of
levels that one of them variables can take, then the ith row ofX(r+1), . . . ,X(c)
will only consist of zeros. Hence, the ith row of L(k) and L˜(k) in Algorithm 1
and the corresponding rows and columns of N(gk) and N˜(gk), k = r+1, . . . , c,
or g = r+1, . . . , c, in Algorithm 2 will be composed of zeros. Since this leads
to dividing zero by zero in (2.2) and (2.5), no χ2-statistic for the ith variable
will be available.
A solution to this problem is to set ˜`
(k)
ig = max
{
1, ˜`
(k)
ig
}
in L˜(k), and
n˜
(k)
ih = max
{
1, n˜
(k)
ih
}
in N˜(k).
3.4 Computation of Similarity Measures
As mentioned in Section 3.1, in the actual realization of Algorithm 2, an
m(m − 1)/2 × c2 matrix M containing all m(m − 1)/2 contingency tables
of the pairwise comparisons is constructed that enables the computation of
similarity measures such as the simple matching coefficient
SM =
1
n
c∑
k=1
nkk.
Since the
(
(k − 1)c + h)th column of M comprises the entry nkh of each
of the contingency tables, the vector rS consisting of SM for any pairwise
comparison of two variables is given by
rS =
Md
M1c2
,
where d is a vector of length c2 with elements
dh =
1, if h ∈
{
a : a = (k − 1)c+ k, k = 1, . . . , c}
0 otherwise
.
9
4 Discussion
In this paper, we have presented approaches based on matrix algebra for the
simultaneous computations of thousands of Pearson’s χ2-statistics.
Table 4.1 shows that the computation time of Algorithm 1 depends only
slightly on the number of classes. However, the number of levels that a
variable can take has a high influence on the computation time time. This
is not very surprising, as the larger c, the more matrices X(k) have to be
constructed and evaluated in the determination of the χ2-statistics.
TABLE 4.1. Computation times of Algorithm 1 for different numbers m of
variables, numbers c of levels a variable can take, and numbers r of classes to
which the n = 200 observations belong.
r = 2, r = 2, r = 2, r = 3, r = 6,
m c = 3 c = 5 c = 10 c = 3 c = 3
100 < 0.01 0.01 0.01 < 0.01 0.01
1,000 0.05 0.07 0.15 0.05 0.05
10,000 0.63 1.03 2.04 0.64 0.62
100,000 6.16 9.98 61.82 6.18 6.46
As Table 4.2 and 4.3 reveal, using Algorithms 1 and 2 lead to a substantial
decreases of the computation time in comparison to one-by-one determina-
tions of the χ2-statistics. Not very surprisingly, the more χ2-statistics, the
higher is the factor by which the computation is accelerated. But even if the
number of variables is small, the algorithms will be about 15 times faster.
Algorithm 1 is used in version 1.10.0 and later of the R package siggenes
available at http://www.bioconductor.org, the web page of the Bioconduc-
tor project (Gentleman et al., 2004), such that the computation time of both
10
TABLE 4.2. Computation times of both Algorithm 1 and the individual calcula-
tion of Pearson’s χ2-statistics for different numbers m of variables and numbers n
of observations. Each variable can take c = 3 levels, and each observation belongs
to one of r = 2 classes.
Algorithm 1 Individual
m n = 200 n = 1, 000 n = 200 n = 1, 000
50 < 0.01 0.01 0.13 0.16
100 < 0.01 0.02 0.26 0.32
1,000 0.05 0.40 2.64 3.35
10,000 0.63 2.39 26.74 34.42
100,000 6.16 – 274.96 –
TABLE 4.3. Computation times of both Algorithm 2 and the individual cal-
culation of Pearson’s χ2-statistic for testing each pair of m variables if they are
independent, where each variable exhibits c = 3 levels, and the number of obser-
vations is n = 1, 000.
m Algorithm 2 Individual
10 0.01 0.15
50 0.07 4.25
100 0.33 17.32
200 1.22 70.06
500 7.79 474.22
11
SAM (Significance Analysis of Microarrays; Tusher et al., 2001) and EBAM
(Empirical Bayes Analysis of Microarrays; Efron et al., 2001) applied to cat-
egorical data (Schwender, 2005, 2007) is reduced. Both Algorithm 1 and 2
are implemented in the R function rowChisqStats contained in the package
scrime that will be available soon at http://www.bioconductor.org and
http://cran.r-project.org.
Acknowledgements
The financial support of the Deutsche Forschungsgemeinschaft (SFB 475,
“Reduction of Complexity in Multivariate Data Structures”) is gratefully
acknowledged.
References
Dudoit, S., Shaffer, J. P., and Boldrick, J.C. (2003). Multiple Hypothesis
Testing in Microarray Experiments. Statistical Science, 18(1), 71–103.
Efron, B., Tibshirani, R., Storey, J.D., and Tusher, V. (2001). Empirical
Bayes Analysis of a Microarray Experiment. Journal of the American Statistical
Association, 96, 1151–1160.
Gentleman, R.C., Carey, V.J., Bates, D.M., Bolstad, B., Dettling, M.,
Dudoit, S., Ellis, B., Gautier, L., Ge, Y., Gentry, J., Hornik, K.,
Hothorn, T., Huber, W., Iacus, S., Irizarry, R., Leisch, F., Li, C.,
Maechler, M., Rossini, A.J., Sawitzki, G., Smith, C., Smyth, G., Tier-
ney, L., Yang, J.Y.H., and Zhang, J. (2004). Bioconductor: Open Software
Development for Computational Biology and Bioinformatics. Genome Biology,
5, R80.
12
Ihaka, R., and Gentleman, R. (1996). R: A Language for Data Analysis and
Graphics. Journal of Computational and Graphical Statistics, 5, 299–314.
Mu¨ller, T., Selinski, S., and Ickstadt, K. (2005). Cluster Analysis: A Com-
parison of Different Similarity Measures for SNP Data. Technical Report, SFB
475, Department of Statistics, University of Dortmund.
Schwender, H. (2005). Modifying Microarray Analysis Methods for Categorical
Data – SAM and PAM for SNPs. In Weihs, C., and Gaul, W. (eds.), Classifica-
tion – The Ubiquitous Challenge. Springer, Heidelberg, 370–377.
Schwender, H. (2007). Empirical Bayes Analysis of Single Nucleotide Polymor-
phisms. In Preparation.
Tusher, V., Tibshirani, R., and Chu, G. (2001). Significance Analysis of Mi-
croarrays Applied to the Ionizing Radiation Response. Proceedings of the Na-
tional Academy of Sciences, 98, 5116–5124.
13


A note on the simultaneous computation of thousands of Pearson’s Chi^2-statistics

https://eldorado.tu-dortmund.de/bitstream/2003/24435/1/tr19-07.pdf

A note on the simultaneous computation of thousands of Pearson’s Chi^2-statistics

Abstract

Similar works

Full text

Available Versions

Eldorado - Ressourcen aus und für Lehre, Studium und Forschung