Schwender, Holger

English

Research Papers in Economics

econstor
www.econstor.eu
Der Open-Access-Publikationsserver der ZBW – Leibniz-Informationszentrum Wirtschaft
The Open Access Publication Server of the ZBW – Leibniz Information Centre for Economics
Nutzungsbedingungen:
Die ZBW räumt Ihnen als Nutzerin/Nutzer das unentgeltliche,
räumlich unbeschränkte und zeitlich auf die Dauer des Schutzrechts
beschränkte einfache Recht ein, das ausgewählte Werk im Rahmen
der unter
→  http://www.econstor.eu/dspace/Nutzungsbedingungen
nachzulesenden vollständigen Nutzungsbedingungen zu
vervielfältigen, mit denen die Nutzerin/der Nutzer sich durch die
erste Nutzung einverstanden erklärt.
Terms of use:
The ZBW grants you, the user, the non-exclusive right to use
the selected work free of charge, territorially unrestricted and
within the time limit of the term of the property rights according
to the terms specified at
→  http://www.econstor.eu/dspace/Nutzungsbedingungen
By the first use of the selected work the user agrees and
declares to comply with these terms of use.
zbw
Leibniz-Informationszentrum Wirtschaft
Leibniz Information Centre for Economics
Schwender, Holger
Working Paper
A note on the simultaneous
computation of thousands of Pearson's
X2-Statistics
Technical Report // Sonderforschungsbereich 475, Komplexitätsreduktion in Multivariaten
Datenstrukturen, Universität Dortmund, No. 2007,19
Provided in cooperation with:
Technische Universität Dortmund
Suggested citation: Schwender, Holger (2007) : A note on the simultaneous computation
of thousands of Pearson's X2-Statistics, Technical Report // Sonderforschungsbereich 475,
Komplexitätsreduktion in Multivariaten Datenstrukturen, Universität Dortmund, No. 2007,19,
http://hdl.handle.net/10419/25004A 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
identiﬁcation of single nucleotide polymorphisms (SNPs) showing a
distribution that diﬀers 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.
11 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 identiﬁca-
tion of SNPs showing a distribution that diﬀers 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 ﬁrst goal is to reduce the
number of SNPs to a better manageable size.
Since the latter goal is similar to the ﬁrst 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 diﬀers between several groups is
computed. The higher this score, the more likely it is that the corresponding
SNP diﬀers 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 X
g=1
c X
k=1
(ngk − ˜ ngk)
2
˜ ngk
, (1.1)
where ngk and ˜ ngk 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 speciﬁes how many of the n =
P
g,k ngk =
P
g,k ˜ ngk
observations in the gth class showing the kth genotype.
A solution to the other problem, i.e. the detection of groups composed of
2related 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 coeﬃcient
PC =
s
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 diﬀerent 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
ﬂexible matching coeﬃcients (M¨ uller 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 diﬀering 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
3are 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 diﬀers 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 = k
0 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
01n, (2.1)
4where 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 diﬀers between the groups speciﬁed 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 X
g=1
c X
k=1
n2
gk
˜ ngk
− 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.
5Algorithm 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 = k
0 otherwise
.
2. For g,k = 1,...,c, compute N(gk) = X(g)X(k)0 and
˜ N
(gk) =
1
n
X
(g)1n1
0
nX
(k)0. (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 X
g=1
c X
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).
73.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
x
A
ij =

 
 
1, xij is not missing
0, xij is a missing value
,
and replace n in (2.3) by the vector
n =

  

n1
. . .
nm

  

= X
A1n.
In (2.1), 1
n is substituted by the m × r matrix
N
den =

  

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 = X
A  
X
A0
.
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)  
XA0 such that (2.4) becomes
˜ N
(gk) =
Z(g) ∗
 
Z(k)0
N
.
83.3 Diﬀerent Numbers of Levels
If the ith variable exhibits r < c levels, where c is the maximum number of
levels that one of the m variables can take, then the ith row of X(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
n
1, ˜ `
(k)
ig
o
in ˜ L(k), and
˜ n
(k)
ih = max
n
1, ˜ n
(k)
ih
o
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 coeﬃcient
SM =
1
n
c X
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
.
94 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 inﬂuence 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 diﬀerent 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
10TABLE 4.2. Computation times of both Algorithm 1 and the individual calcula-
tion of Pearson’s χ2-statistics for diﬀerent 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
11SAM (Signiﬁcance 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 ﬁnancial 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.
12Ihaka, R., and Gentleman, R. (1996). R: A Language for Data Analysis and
Graphics. Journal of Computational and Graphical Statistics, 5, 299–314.
M¨ uller, T., Selinski, S., and Ickstadt, K. (2005). Cluster Analysis: A Com-
parison of Diﬀerent 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.), Classiﬁca-
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). Signiﬁcance Analysis of Mi-
croarrays Applied to the Ionizing Radiation Response. Proceedings of the Na-
tional Academy of Sciences, 98, 5116–5124.
13

Bioconductor: Open Software Development for Computational Biology and Bioinformatics.

Cluster Analysis: A Comparison of Diﬀerent Similarity Measures for SNP Data.

Empirical Bayes Analysis of a Microarray Experiment.

Empirical Bayes Analysis of Single Nucleotide Polymorphisms. In Preparation.

Modifying Microarray Analysis Methods for Categorical Data

Multiple Hypothesis Testing in Microarray Experiments.

R: A Language for Data Analysis and Graphics.

Signiﬁcance Analysis of Microarrays Applied to the Ionizing Radiation Response.

A note on the simultaneous computation of thousands of Pearson's X2-Statistics

http://econstor.eu/bitstream/10419/25004/1/550475796.PDF

A note on the simultaneous computation of thousands of Pearson's X2-Statistics

Abstract

Similar works

Full text

Available Versions

Research Papers in Economics