BMC Bioinformatics, vol.10, no.1, pp.S52, 2009.

1

Mutual Information Estimation Reveals Global
Associations between Stimuli and Biological Processes
Taiji Suzuki (s-taiji@stat.t.u-tokyo.ac.jp)
Department of Mathematical Informatics, The University of Tokyo
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
Masashi Sugiyama (sugi@cs.titech.ac.jp)
Department of Computer Science, Tokyo Institute of Technology
2-12-1, O-okayama, Meguro-ku, Tokyo, 152-8552, Japan
Takafumi Kanamori (kanamori@is.nagoya-u.ac.jp)
Department of Computer Science and Mathematical Informatics,
Nagoya University
Furocho, Chikusaku, Nagoya 464-8603, JapanJun Sese (sesejun@is.ocha.ac.jp)
Department of Information Science, Ochanomizu University
2-1-1 Ohtsuka, Bunkyo-ku, Tokyo 112-8610, Japan
Abstract
Background: Although microarray gene expression analysis has become popular,
it remains diﬃcult to interpret the biological changes caused by stimuli or variationof conditions. Clustering of genes and associating each group with biological functions are often used methods. However, such methods only detect partial changes
within cell processes. Herein, we propose a method for discovering global changes
within a cell by associating observed conditions of gene expression with gene functions.
Results: To elucidate the association, we introduce a novel feature selection method
called Least-Squares Mutual Information (LSMI), which computes mutual information without density estimaion, and therefore LSMI can detect nonlinear associations within a cell. We demonstrate the eﬀectiveness of LSMI through comparison
with existing methods. The results of the application to yeast microarray datasets
reveal that non-natural stimuli aﬀect various biological processes, whereas others
are no signiﬁcant relation to speciﬁc cell processes. Furthermore, we discover that
biological processes can be categorized into four types according to the responses ofvarious stimuli: DNA/RNA metabolism, gene expression, protein metabolism, andprotein localization.
Conclusions: We proposed a novel feature selection method called LSMI, and applied LSMI to mining the association between conditions of yeast and biologicalprocesses through microarray datasets. In fact, LSMI allows us to elucidate theglobal organization of cellular process control.

Mutual Information Estimation Reveals Global Associations

2

Background
Advances in microarray technologies enable us to explore the comprehensive dynamics
of transcription within a cell. The current problem is to extract useful information froma massive dataset. The primarily used approach is clustering. Cluster analysis reveals
variations of gene expression and reduces the complexity of large datasets. However,additional methods are necessary to associate genes in each cluster with genetic function
using GO term ﬁnder [3], or to understand stimuli related to speciﬁc cellular status.
However, these clustering-association strategies cannot detect global cell statuschanges because of the division of clusters. Some stimuli activate a speciﬁc pathway,
although others might change overall cellular processes. Understanding the eﬀect of stimuli in cellular processes directly, in this paper, we introduce a novel feature selection
method called Least-Squares Mutual Information (LSMI), which selects features using
mutual information without density estimation. Mutual information has been utilized tomeasure distances between gene expressions [15]. To compute the mutual information inexisting methods, density estimation or discritization is required. However, the estimationof gene expression is diﬃcult because we have little knowledge about density function ofgene expression proﬁle. LSMI oﬀers an analytic-form solution and avoid the estimation.
Feature selection techniques are often used in gene expression analysis [16]. Actually,
LSMI has three advantages compared to existing methods: capability of avoiding density
estimation which is known to be a hard problem [18], availability of model selection, andfreedom from a strong model assumption. To evaluate the reliability of ranked featuresusing LSMI, we compare receiver operating characteristic (ROC) curves [14] to those ofexisting methods: kernel density estimation (KDE) [20, 6], k-nearest neighbor (KNN) [13],Edgeworth expansion (EDGE) [11], and Pearson correlation coeﬃcient (PCC). Thereby,
we certify that our method has better performance than the existing methods in predictionof gene functions about biological processes. This fact implies that features selected using
our method reﬂect biological processes.
Using the ranked features, we illustrate the associations between stimuli and biologicalprocesses according to gene expressions. Results show that stimuli damage essential processes within a cell, causing association with some cellular processes. From the responseto stimuli, biological processes are divisible into four categories: DNA/RNA metabolic
processes, gene expression, protein metabolic processes, and protein localization.

Results
Approach—Mutual Information Detection
In this study, we detect underlying dependencies between gene expressions obtained bygroups of stimuli and gene functions. The dependencies are studied in various machinelearning problems such as feature selection [8, 22] and independent component analysis [4].
Although classical correlation analysis would be useful for these problems, it cannot detect
nonlinear dependencies with no correlation. On the other hand, mutual information (MI),Mutual Information Estimation Reveals Global Associations

3

which plays an important role in information theory [5], enables us to detect general
nonlinear dependencies. Let x and y be a set of gene expressions and a set of known genefunctions. A variant of MI based on the squared loss is deﬁned by)2∫∫ (pxy (x, y)− 1 px (x)py (y)dxdy.
(1)
Is (X, Y ) :=px (x)py (y)Note that Is vanishes if and only if x and y are independent. The use of MI allows us todetect no correlation stimulus with a speciﬁc gene function or process.
Estimating MI is known to be a diﬃcult problem in practice [22, 13, 11]. Herein, we
propose LSMI, which does not involve density estimation but directly models the density
ratio:
pxy (x, y)w(x, y) :=.
px (x)py (y)Given a density ratio estimator w(x,b y), squared loss MI can be simply estimated byn1 ∑bIs (X, Y ) = 2(w(xb i , y j ) − 1)2 .
n i,j=1

Mathematical deﬁnitions related to LSMI are provided in the Methods section. LSMIoﬀers an analytic-form solution, which allows us to estimate MI in a computationallyvery eﬃciently manner.
It is noteworthy that x includes a multi-dimensional vector. In fact, LSMI can handle agroup of stimuli, although generic correlation indices such as Pearson correlation between
parameters and target value are calculated independently. Therefore, we can elucidate
which type of stimulus has no dependency to biological processes using LSMI.

Datasets and Feature Selection
In this section, we ﬁrst prepare datasets to show the association between stimuli andbiological process, and introduce feature selection using the datasets.
Biological Process
We compute mutual information between gene expression values grouped by stimuli andclass of genes’ biological processes. As the class, we use biological process terms in Gene
Ontology (GO) categorization [1]. We select GO terms associated with more than 800
and less than 2,000 genes because terms having a small number of genes only describe afraction of the cell status, whereas terms having a large number of genes indicate functions
associated with almost all genes in yeast. Actually, GO has a directed acyclic graph(DAG) structure, and each term has child terms. The GO terms are classiﬁed into three
categories; we use only biological process terms to identify the changes within a cell. Using
this method, we select 12 GO terms.

Mutual Information Estimation Reveals Global Associations

4

Gene Expression Profiles
The gene expression proﬁle is the best comprehensive dataset to associate stimuli andbiological processes. We use two diﬀerent microarray datasets. One is of 173 microarraydata under stress conditions of various types [7]. We categorize the 173 stress conditions
into 29 groups based on the type of condition such as heat shock, oxidizing condition, etc.
The other is of 300 microarray data under gene-mutated conditions [10]. We categorize
the genes into 146 groups based on associated GO terms. We use only the GO terms
which are associated with 1,500 genes or fewer. We also use child terms on a GO layered
structure if the term has more than 1,200 genes. When one gene belongs to multiple GOterms, we classify the gene into the the classiﬁcation whose number of associated genes issmallest.
In both proﬁles, we remove genes whose expression values are obtained from fewer
than 30% of all observed conditions. All missing values are ﬁlled out by the average of all
the expression values.
Feature Selection using LSMI
We use a novel feature selection method called LSMI, which is based on MI, to associate stimuli with cellular processes. Here we consider the forward feature-group addition
strategy, i.e., a feature-group score between each input feature-group and output cellular
process is computed. The top m feature-groups are used for training a classiﬁer. We predict 12 GO terms independently. We randomly choose 500 genes from among 6, 116 genes
on the stress condition dataset for feature-group selection and for training a classiﬁer; therest are used for evaluating the generalization performance. For using the gene-mutatedexpression dataset, we select 500 genes from among 6, 210 genes. We repeat this trial 10
times. For classiﬁcation, we use a Gaussian kernel support vector machine (GK-SVM)[18], where the kernel width is set at the median distance among all samples and theregularization parameter is ﬁxed at C = 10. We explain the eﬃciency of feature selectionof LSMI in the Discussion section.

Results
The association between stress conditions and biological processes in GO terms is shownin Fig.1. Each row and column respectively indicate a group of conditions and a GO term.
Row and column dendrograms are clustering results by the Ward method according tocell values. Each cell contains an average ranking over 10 trials by LSMI. The red cell
denotes that the parameter has a higher rank; that is, the parameter has association withthe target GO term. A blue cell denotes that the parameter has a lower rank.
As shown in this ﬁgure, conditions are divided into two groups. Almost all conditionsin the upper cluster have higher rank, whereas those in a lower cluster have higher rank
only under speciﬁc conditions. The conditions in the upper cluster include strong heatshocks, dithiothreitol (DTT) exposure, nitrogen depletion, and diamide treatments, whichMutual Information Estimation Reveals Global Associations

5

Over expression + 37C{17,21,..,33}C to 37C
DTT <= 3h
Nitrogen depletion
Amino acid starvations

37C to 25C
DiamideStationary phase 30C (<= 12h)
DTT > 3h
Stationary phase 25C (<= 12h)
Stationary phase 30C (> 12h)
Diauxic shift
Stationary phase 25C (> 12d)Over expression
29C +1m sorbitol to 33C

(B)

G
O
:0
0

44

26

0:

ce
l

G
O
:0
0
lu 195lar_ 3 8 :mGa c p ro
O
te
ro
:0
in
m
00
_
ol
90
ec me
59
ta
ul
G
bo
:m
e_
O
lic
m
:0
ac
e
0
ro
ta _pr
10
m
o
bo
4
ol
67
lic ces
ec
G
:g
_p
s
O
u
e
r
l
:0
e_
oc
ne
00
bi
es
_
90
ex
os
G
s
O
pr
:0 58: ynt
es
he
06
bi
si
o
ti
50
G
on
O
07 syn c_p:0: b the ro
0
G
ce
51
G
O
i
t
O
:0
s
23
G olog ic_
00 :00
O
pr s
4:
ic
:
o
69
43
0
al
e
c
0
s
_
e
G
51
ta
r
s
O 96: 283
bl
17 egu s:0is
9:0 0 o rg : b i
la
hm
o
a
t
lo
61
p
n
en
ca ion39 elle oly
liz:n_o me t_o
at
f
_a uc
_l
r
r
io
oc
n
nd le gan _m
et
al
iz
_n ob
a
i
a
z
a
b
at
tio
G uc seoiol
O le ,_
n:0 ic n n_a ic_
01 _a ucnd pro60 cid lec_esbi
70 _m os
: R e ide oge s
N tab ,_n ne
A_ o usm lic_ cle iset p ota b ro i d
ol ce e
ic ss
_p
ro
ce
ss

29C to 33C
Over expression + H2O2
Sorbitol
Ethanol
Menadione
H2O2
Monosaccharides
Polysaccharides(A)

Figure 1: Matrix of stress conditions (rows) versus biological processes (columns). Red
cells have higher correlation.
are non-natural conditions. The result reveals that non-natural conditions change overallcellular processes.
The GO term clusters are divided into three groups: DNA/RNA metabolism (right),
localization of protein (middle), and others (left). The leftmost cluster contains bio synthesis, gene expression process, and protein metabolic process. From this ﬁgure, nucleic
acid metabolism processes are inferred to be independent from amino acid metabolism
processes. We will conﬁrm the independence and consider the division of clusters by using
other dataset later.
We herein investigate the details of diﬀerence among DNA metabolic process, protein
metabolic process and localization of proteins. Under an overexpression condition indicated by sign (A) in Fig.1, DNA/RNA metabolisms show no correlation with expressions

GO:0016070, RNA_metabolic_process

GO:0044260, cellular_macromolecule_metabolic_process

GO:0019538, protein_metabolic_process

GO:0009058, biosynthetic_process

GO:0009059, macromolecule_biosynthetic_process

GO:0065007, biological_regulation

GO:0051234, establishment_of_localization

GO:0051179, localization

GO:0006996, organelle_organization_and_biogenesis

GO:0010467, gene_expression

GO:0043283, biopolymer_metabolic_process

GO:0006139, nucleobase,_nucleoside,_nucleotide_and_nucleic_acid_metabolic_process

Mutated gene sets

Mutual Information Estimation Reveals Global Associations
6

(III)

(II)

(I)

Figure 2: Overview: a matrix of mutated gene groups (rows) versus biological processes(columns)

Mutual Information Estimation Reveals Global Associations(I) High correlations with any biological process

7

Mutated Gene Groups
chromosome organization and biogenesis
cellular homeostasis
cellular component assembly
cellular localization
cell cycle processtranslation
lipid biosynthetic processexternal encapsulating structure organization
and biogenesis
cell development
HU
growthmitochondrion organization and biogenesis
RNA processing
chromosome segregation
glycoprotein metabolic process(II) Association with DNA/RNA metabolic process

peptide metabolic process
Cycloheximide
yel033w
arg5,6(III) Signiﬁcantly no association with DNA/RNA metabolic process
yer084w
yel010w
yer041w
yhl045w
yor009w
ymr141c

GO:0044260, cellular_macromolecule_metabolic_process

GO:0019538, protein_metabolic_process

GO:0009058, biosynthetic_process

GO:0009059, macromolecule_biosynthetic_process

GO:0065007, biological_regulation

GO:0051234, establishment_of_localization

GO:0051179, localization

GO:0006996, organelle_organization_and_biogenesis

GO:0010467, gene_expression

GO:0043283, biopolymer_metabolic_process

GO:0016070, RNA_metabolic_process

GO:0006139, nucleobase,_nucleoside,_nucleotide_and_nucleic_acid_metabolic_process

ymr147w
yml003w

Figure 3: Sub-matrices of the full mapof genes belonging to over-expression genes. This ﬁnding of no correlation is one advantage of LSMI. The menadione (vitamin K) exposure condition indicated by (B) in Fig.1 isassociated with localization of proteins. Menadione supplementation causes high toxicity;such toxicity might result from the violation of protein localizations.
Next, we compute the association using expressions of gene mutants. The results areshown in Fig.2. The stimulus can be categorized into two parts: high association underalmost all processes and under particular conditions. The division is the same because

Mutual Information Estimation Reveals Global Associations

8of stress condition associations. The GO terms also categorize three parts: DNA/RNA
metabolic processes, protein metabolic processes, and localization. In this experiment,
GO terms “gene expression” (GO:0010467) and “organelle organization and biogenesis”(GO:0006996) are in the DNA/RNA metabolic process cluster, although they are classiﬁed in protein metabolic processes cluster under stress conditions in Fig.1. Because
the both divisions are close to ancestor division, we can conclude that the cluster about
gene expression exists. From these results, GO terms are divisible into four categories:
DNA/RNA metabolic process, protein metabolic process, localization, and gene expression.
In Fig.3, we present details of three clusters in Fig.2. In fact, Fig.3(I) presents a cluster
whose members are correlated with any biological process. Furthermore, the functions of
the mutated genes are essential processes for living cells, such as cellular localization, cell
cycle, and growth. This result might indicate that the upper half stimulus in Fig.1 destroys
the functions of these essential genes. Furthermore, Fig.3(II) includes the groups of genes
associated with DNA/RNA metabolic processes. In this cluster, YEL033W/MTC1 is agene with unknown function and is predicted to have a metabolic role using protein–protein interaction [17]. Our clustering result indicates that YEL033W would have somerelation with metabolism, especially methylation (methylation is an important part ofthe one-carbon compound metabolic process). We show genes which have no signiﬁcantassociation with DNA/RNA metabolic processes in Fig.3(III). In the cluster, all genes
except AQY2 are of unknown function. No correlation clusters cannot be found by existing
methods. Our result might provide clues to elucidate these genes’ functions.

Discussion
A common analytical ﬂow of the expression data is ﬁrst clustering and then associatingclusters with GO terms or pathways. Although clustering reduces the complexity of large
datasets, the strategy might fail to detect changes of entire genes within a cell such as
metabolic processes.
To interpret such gene expression changes, gene set enrichment analysis [21] has been
proposed. This method treats microarrays independently. Therefore, housekeeping genes
are often ranked highly. When gene expressions under various conditions are available, ourmethod would show us the better changes of cellular processes because of the comparison
between groups of conditions. The module map [19] gives a global association between
a set of genes and a set of conditions. However, this method requires important changesof gene expressions because it uses hypergeometric distributions to compute correlations.
Our correlation index is based on MI. Therefore, we can detect nonlinear dependencies
with no correlation. An example is depicted in Fig.3(III).
The characteristics of LSMI and existing MI estimators are presented in Table 1.
Detail comparisons are described in the Methods section. The kernel density estimator(KDE) [20, 6] is distribution-free. Model selection is possible by likelihood cross-validation(LCV). However, a hard task of density estimation is involved. Estimation of the entropies

Mutual Information Estimation Reveals Global Associations

9

Table 1: Relation between existing and proposed MI estimators. If the order of the
Edgeworth expansion is regarded as a tuning parameter, model selection of EDGE is
expected to be ‘Not available’.
KDE
KNN
EDGE
LSMI

Density estimation
Involved
Not involved
Not involved
Not involved

Model selection
Available
Not available
Not necessary
Available

Distribution
Free
Free
Nearly normalFreeusing k-nearest neighbor (KNN) samples [13] is distribution-free and does not involve density estimation directly. However, no model selection method exists for determining thenumber of nearest neighbors. Edgeworth expansion (EDGE) [11] does not involve density estimation or any tuning parameters. However, it is based on the assumption thatthe target distribution is close to the normal distribution. On the other hand, LSMI
is distribution-free; it involves no density estimation, and model selection is possible bycross-validation (CV). Therefore, LSMI overcomes limitations of the existing approaches.
Within a cell, most processes have a nonlinear relation such as enzyme eﬀects and feedback loops. The lack of one advantage might cause diﬃculty of application to biologicaldatasets. By virtue of these advantages, LSMI can detect correlation or independence
between features of complex cellular processes.
To investigate the eﬃciency of feature selection, we compare areas under the curve(AUCs) with LSMI (CV), KDE(LCV), KNN(k) for k = 1, 5, EDGE, and PCC. Detailsof these methods are described in the Methods section. Figure 4 depicts AUCs for 12
GO term classiﬁcations. The x-axis shows the number of stimulus groups used for theprediction. The y-axis means averaged AUC over 10 trials, where AUCs are calculated
as the area under the receiver operating characteristic (ROC) curve, which is often used
for diagnostic tests. Each ﬁgure shows AUC curves calculated using the six methods.
In the AUC ﬁgures, the higher curves represent better predictions. For example,
Fig.4(a) shows that LSMI is the highest position, which means that LSMI achieves thebest performance among the six methods. In Figs.4(b) and 4(d), KNN(1) and KNN(5),which are denoted by the light blue and dotted light blue lines, have the best performance.
However, in Figs.4(i), 4(j) and 4(l), averaged AUCs of KNN using numerous groups
are high, whereas the AUCs using small and few groups are low. No systematic model
selection strategies exist for KNN and therefore KNN would be unreliable in practice.
Figure 4(c) depicts that EDGE, which is indicated by the light green line, has the highest
AUC. In fact, EDGE presumes the normal distribution. Consequently, it works well only
on a few datasets. From these ﬁgures, LSMI indicated by the blue line appears to be thebest feature selection method.

Mutual Information Estimation Reveals Global Associations

LSMI

KDE

KNN(1)

KNN(5)

10

PCC

EDGE

0.64

0.64
0.62
0.6
0.58
0.56
1

0.64

Averaged AUC over 10 trials

Averaged AUC over 10 trials

Averaged AUC over 10 trials

0.66

0.62
0.6
0.58

4

Number of Groups

10

29

1

(a) nucleotide and nucleic
acid metabolic process2

4

Number of Groups

1029(b) organelle organization
and biogenesis

2

4

Number of Groups

1029(c) biosynthetic process

0.64
0.62
0.6
0.58
0.56

Averaged AUC over 10 trials

0.7
Averaged AUC over 10 trials

Averaged AUC over 10 trials

0.58

1

0.66

1

0.6

0.56

0.56
2

0.62

0.7

0.65

0.65

0.6

0.6
2

4

Number of Groups

10

29

1

(d) macromolecule biosynthetic process2

4

Number of Groups

10

0.55
1

29

(e) gene expression

2

4

Number of Groups

1029(f) RNA metabolic process

0.66
0.64
0.62
0.6
0.58
4

Number of Groups

10

29

0.58
0.56

1

0.57

Averaged AUC over 10 trials

Averaged AUC over 10 trials

0.6

0.54
2

(g) protein metabolic process
0.56
0.55
0.54
0.53
0.52
1

0.62

2

4

Number of Groups

10(j) localization

29

0.66
0.64
0.62
0.6
0.58
0.56

2

4

Number of Groups

10

29

1

2

4

Number of Groups

1029(h) biopolymer metabolic
process(i) cellular macromolecule
metabolic process

0.58

0.58
Averaged AUC over 10 trials

1

0.64

Averaged AUC over 10 trials

Averaged AUC over 10 trials

Averaged AUC over 10 trials

0.68

0.56
0.54
0.52

1

2

4

Number of Groups

1029(k) establishment of localization

0.57
0.56
0.55
0.54
0.53
0.52
1

2

4

Number of Groups

1029(l) biological regulation

Figure 4: Classiﬁcation error against the number of feature groups for the yeast cell
datasets.

Mutual Information Estimation Reveals Global Associations

11

Conclusions
We provided a global view of the associations between stimuli and changes of biologicalprocesses based on gene expression proﬁles. The association is generally diﬃcult to use
for making models because of nonlinear correlation. To cope with this problem, we introduced a novel feature selection method called LSMI, which uses MI and can be computed
eﬃciently. In comparison to other feature selection methods, LSMI showed better AUCs
in prediction of biological process functions. Consequently, our feature selection results
would be more reliable than those obtained using the other methods. We calculated
the association between stimuli and GO biological process terms using gene expression
proﬁles. The result revealed that the stimuli are categorized into four types: related to
DNA/RNA metabolic process, gene expression, protein metabolic process, and protein
localization. LSMI enabled us to reveal the global regulation of cellular processes fromcomprehensive transcription datasets.

Methods
Mutual Information Estimation
A naive approach to estimating MI is to use a KDE [20, 6], i.e., the densities pxy (x, y),px (x), and py (y) are separately estimated from samples and the estimated densities areused for computing MI. The band-width of the kernel functions could be optimized basedon likelihood cross-validation (LCV) [9], so there remains no open tuning parameter inthis approach. However, density estimation is known to be a hard problem [18] andtherefore the KDE-based method may not be so eﬀective in practice.
An alternative method involves estimation of entropies using KNN. The KNN-based
approach was shown to perform better than KDE [12], given that the number k is chosen
appropriately—a small (large) k yields an estimator with small (large) bias and large(small) variance. However, appropriately determining the value of k is not straightforwardin the context of MI estimation.
Here, we propose a new MI estimator that can overcome the limitations of the existing
approaches. Our method, which we call Least-Squares Mutual Information (LSMI), does
not involve density estimation and directly models the density ratio:
w(x, y) :=pxy (x, y).
px (x)py (y)(2)

The solution of LSMI can be computed by simply solving a system of linear equations.
Therefore, LSMI is computationally very eﬃcient. Furthermore, a variant of crossvalidation (CV) is available for model selection, so the values of tuning parameters suchas the regularization parameter and the kernel width can be adaptively determined in an
objective manner.

Mutual Information Estimation Reveals Global Associations

12

A New MI Estimator
In this section, we formulate the MI inference problem as density ratio estimation andpropose a new method of estimating the density ratio.
MI Inference via Density Ratio Estimation
Let DX (⊂ Rdx ) and DY (⊂ Rdy ) be the data domains and suppose we are given n
independent and identically distributed (i.i.d.) paired samples{(xi , y i ) | xi ∈ DX , y i ∈ DY }ni=1drawn from a joint distribution with density pxy (x, y). Let us denote the marginal densities of xi and y i by px (x) and py (y), respectively. The goal is to estimate squared-loss
MI deﬁned by Eq.(1).
Our key constraint is that we want to avoid density estimation when estimating MI.
To this end, we estimate the density ratio w(x, y) deﬁned by Eq.(2). Given a density
ratio estimator w(x,b y), MI can be simply estimated byn1 ∑bIs (X, Y ) = 2(w(xb i , y j ) − 1)2 .
n i,j=1

We model the density ratio function w(x, y) by the following linear model:wbα (x, y) := α⊤ φ(x, y),where α = (α1 , α2 , . . . , αb )⊤ are parameters to be learned from samples, ⊤ denotes thetranspose of a matrix or a vector, andφ(x, y) = (φ1 (x, y), φ2 (x, y), . . . , φb (x, y))⊤are basis functions such thatφ(x, y) ≥ 0b for all (x, y) ∈ DX × DY .
0b denotes the b-dimensional vector with all zeros. Note that φ(x, y) could be dependent
on the samples {xi , y i }ni=1 , i.e., kernel models are also allowed. We explain how the basis
functions φ(x, y) are chosen in the later section.
A Least-squares Approach to Direct Density Ratio Estimation
We determine the parameter α in the model w
bα (x, y) so that the following squared error
J0 is minimized:∫∫
1(wbα (x, y)−w(x, y))2px (x)py (y)dxdy
J0 (α) :=
2
∫∫
∫∫
1
2
=
w
bα (x, y) px (x)py (y)dxdy −wbα (x, y)pxy (x, y)dxdy + C,2

Mutual Information Estimation Reveals Global Associations

13∫∫
where C = 12
w(x, y)pxy (x, y)dxdy is a constant and therefore can be safely ignored.
Let us denote the ﬁrst two terms by J:
1
J(α) := J0 (α) − C = α⊤ Hα − h⊤ α,
2
where

∫∫
H :=
∫∫
h :=

φ(x, y)φ(x, y)⊤ px (x)py (y)dxdy,
φ(x, y)pxy (x, y)dxdy.

Approximating the expectations in H and h by empirical averages, we obtain thefollowing optimization problem:
[
]
⊤
1 ⊤c
⊤
b
e := argminαα Hα − h α + λα α ,
(3)
2
α∈Rb
where we included a regularization term λα⊤ α andn∑c := 1H
φ(xi , y j )φ(xi , y j )⊤ ,n2 i,j=1n∑b := 1hφ(xi , y i ).
n i=1

Diﬀerentiating the objective function (3) with respect to α and equating it to zero, we
can obtain an analytic-form solution:bc + λI b )−1 h,e = (Hαwhere I b is the b-dimensional identity matrix.
We call the above method Least-Squares Mutual Information (LSMI) . Thanks to theanalytic-form solution, the LSMI solution can be computed very eﬃciently.
Convergence Bound
Here, we show a non-parametric convergence rate of the solution of the optimization
problem (3).
Let G be a general set of functions on DX × DY . For a function g (∈ G), let us consider
a non-negative function R(g) such thatsup[g(x, y)] ≤ R(g).
x,yThen the problem (3) can be generalized as[]nn
1 ∑ 2 1∑
gi,j −w
b := argmin
gi,i + λn R(g)2 ,
2
2n
n
g∈G
i,j=1
i=1

Mutual Information Estimation Reveals Global Associations

14

where gi,j := g(xi , y j ). We assume that the true density ratio function w(x, y) is containedin model G and satisﬁes
w(x, y) < M0 for all (x, y) ∈ DX × DY .
We also assume that there exists γ (0 < γ < 2) such that
H[] (GM , ϵ, L2 (pX pY )) = O((M/ϵ)γ ),where
GM := {g ∈ G | R(g) ≤ M }and H[] is the bracketing entropy of GM with respect to the L2 (px py )-norm [25, 24]. Thismeans the function class G is not too much complex.
Then we have the following theorem. Its proof is omitted due to lack of space.
2/(2+γ)Theorem 1 Under the above setting, if λn → 0 and λ−1) thenn = o(n∥w
b − w∥2 = Op (λ1/2n ),where ∥ · ∥2 means the L2 (px py )-norm and Op denotes the asymptotic order in probability.
This theorem is closely related to [23, 2]. The paper [23] considers least squares estimatorsfor nonparametric regression, and related topics can be found in Section 10 of [24].
CV for Model Selection and Basis Function Design
The performance of LSMI depends on the choice of the model, i.e., the basis functions
φ(x, y) and the regularization parameter λ. Here we show that model selection can be
carried out based on a variant of CV.
First, the samples {z i | z i = (xi , y i )}ni=1 are divided into K disjoint subsets {Zk }K
k=1 .
Then a density ratio estimator w
bk (x, y) is obtained using {Zj }j̸=k and the cost J isapproximated using the held-out samples Zk as(K−CV)Jbk=∑w
bk (x′ , y ′ )2
x′ ,y ′ ∈Zk

2n2k

−

∑ w
bk (x′ , y ′ )(x′ ,y ′ )∈Zk

nk

,

∑
where nk is the number of pairs in the set Zk .
is the summation over all∑ x′ ,y′ ∈Zk′′2
combinations of x and y (i.e., nk terms), while (x′ ,y′ )∈Zk is the summation over all
pairs (x′ , y ′ ) (i.e., nk terms). This procedure is repeated for k = 1, 2, . . . , K and its
average Jb(K−CV) is used as an estimate of J:K
1 ∑ b(K−CV)(K−CV)bJ=.
J
K k=1 r

Mutual Information Estimation Reveals Global Associations

15

We can show that Jb(K−CV) is an almost unbiased estimate of the true cost J, wherethe ‘almost’-ness comes from the fact that the number of samples is reduced in the CV
procedure due to data splitting [18].
A good model may be chosen by CV, given that a family of promising model candidates
is prepared. As model candidates, we propose using a Gaussian kernel model:()∥x − uℓ ∥2
φℓ (x, y) = exp −δ(y = v ℓ ),
2σ 2
where
{(uℓ , v ℓ )}bℓ=1
are ‘center’ points randomly chosen from{(xi , y i )}ni=1 .
δ(y = v ℓ ) is a indicator function, which is 1 if y = v ℓ and 0 otherwise.
In the experiments, we ﬁx the number of basis functions at
b = min(100, n),and choose the Gaussian width σ and the regularization parameter λ by CV with gridsearch.

Relation to Existing Methods
In this section, we discuss the characteristics of existing and proposed approaches.
Kernel Density Estimator (KDE)
KDE [20, 6] is a non-parametric technique to estimate a probability density function p(x)from its i.i.d. samples {xi }ni=1 . For the Gaussian kernel, KDE is expressed as
(
)
n
∑
1
∥x − xi ∥2pb(x) =exp −.
n(2πσ 2 )d/2 i=1
2σ 2
The performance of KDE depends on the choice of the kernel width σ and it can be
optimized by likelihood CV as follows [9]: First, divide the samples {xi }ni=1 into K disjoint
subsets {Xk }K
bXk (x) from {Xj }j̸=k and compute its
k=1 . Then obtain a density estimate p
hold-out log-likelihood for Xk :1 ∑log pbXk (x).
|Xk | x∈X
kThis procedure is repeated for k = 1, 2, . . . , K and choose the value of σ such that the
average of the hold-out log-likelihood over all k is maximized. Note that the averageMutual Information Estimation Reveals Global Associations

16

hold-out log-likelihood is an almost unbiased estimate of the Kullback-Leibler divergencefrom p(x) to pb(x), up to an irrelevant constant.
Based on KDE, MI can be approximated by separately estimating the densities
pxy (x, y), px (x) and py (y) using {xi , y i }ni=1 . However, density estimation is known to
be a hard problem and therefore the KDE-based approach may not be so eﬀective inpractice.
k-nearest Neighbor Method (KNN)
Let Nk (i) be the set of k-nearest neighbor samples of (xi , y i ), and let
ϵx (i) := max{∥xi − xi′ ∥ | (xi′ , y i′ ) ∈ Nk (i)},ϵy (i) := max{∥y i − y i′ ∥ | (xi′ , y i′ ) ∈ Nk (i)},nx (i) := #{z i′ | ∥xi − xi′ ∥ ≤ ϵx (i)},ny (i) := #{z i′ | ∥y i − y i′ ∥ ≤ ϵy (i)}.
Then the KNN-based MI estimator is given as follows [13]:n1 1∑bI(X, Y ) = ψ(k) + ψ(n) − −[ψ(nx (i)) + ψ(ny (i))] ,k n i=1

where ψ is the digamma function.
A practical drawback of the KNN-based approach is that the estimation accuracy
depends on the value of k and there seems no systematic strategy to choose the value ofk appropriately.
Edgeworth Expansion (EDGE)
MI can be expressed in terms of the entropies as
I(X, Y ) = H(X) + H(Y ) − H(X, Y ),where H(X) denotes the entropy of X:∫
H(X) := −px (x) log px (x)dx.

Thus MI can be approximated if the entropies above are estimated.
In the paper [11], an entropy approximation method based on the Edgeworth expansion
is proposed, where the entropy of a distribution is approximated by that of the normal
distribution and some additional higher-order correction terms. More speciﬁcally, for ad-dimensional distribution, the entropy is approximated by
H ≈ Hnormal −

d
d
d
∑
1 ∑ 2
1
1 ∑ 2
κi,i,i −
κi,i,j −κ2 ,12 i=1
4 i,j=1,i̸=j72 i,j,k=1,i<j<k i,j,kMutual Information Estimation Reveals Global Associations

17

where Hnormal is the entropy of the normal distribution with covariance matrix equal to
the target distribution and κi,j,k (1 ≤ i, j, k ≤ d) is the standardized third cumulant of
the target distribution. In practice, all the cumulants are estimated from samples.
If the underlying distribution is close to the normal distribution, the above approximation is quite accurate and the EDGE method works very well. However, if the distribution
is far from the normal distribution, the approximation error gets large and therefore the
EDGE method may be unreliable.
In principle, it is possible to include the fourth and even higher cumulants for further
reducing the estimation bias. However, this in turn increases the estimation variance; theexpansion up to the third cumulants would be reasonable.

Competing interests
The authors declare that they have no competing interests.

Authors contributions
TS developed the method, implemented the algorithm and wrote the manuscript. MS
and TK discussed the method and revised the manuscript. JS discussed the method,interpreted the results and wrote the manuscript.

Acknowledgements
This work was partially supported by KAKENHI (Grant-in-Aid for Scientiﬁc Research)on Priority Areas ”Systems Genomics” from the Ministry of Education, Culture, Sports,
Science and Technology of Japan.
T.S. was supported by the JSPS Research Fellowships for Young Scientists.

References[1] M. Ashburner, C. A. Ball, J. A. Blake, D. Botstein, H. Butler, et al. Gene ontology:
tool for the uniﬁcation of biology. the gene ontology consortium. Nature Genetics,
25(1):25–29, May 2000.
[2] L. Birgé and P. Massart. Rates of convergence for minimum contrast estimators.
Probability Theory and Related Fields, 97:113–150, 1993.
[3] E. I. Boyle, S. Weng, J. Gollub, H. Jin, D. Botstein, J. M. Cherry, and G. Sherlock.
Go::termﬁnder—open source software for accessing gene ontology information andﬁnding signiﬁcantly enriched gene ontology terms associated with a list of genes.
Bioinformatics, 20:3710–3715, 2004.

Mutual Information Estimation Reveals Global Associations[4] P. Comon. Independent component analysis, a new concept?
36(3):287–314, 1994.

18

Signal Processing,[5] T. M. Cover and J. A. Thomas. Elements of Information Theory. John Wiley &Sons, Inc., N. Y., 1991.
[6] A. M. Fraser and H. L. Swinney. Independent coordinates for strange attractors frommutual information. Physical Review A, 33(2):1134–1140, 1986.
[7] A. P. Gasch, P. T. Spellman, C. M. Kao, O. Carmel-Harel, M. B. Eisen, G. Storz,
D. Botstein, and P. O. Brown. Genomic expression programs in the response of yeast
cells to environmental changes. Molecular Biology of the Cell, 11(12):4241–4257,
2000.
[8] I. Guyon and A. Elisseeﬀ. An introduction to variable feature selection. Journal of
Machine Learning Research, 3:1157–1182, 2003.
[9] W. Härdle, M. Müller, S. Sperlich, and A. Werwatz. Nonparametric and Semiparametric Models. Springer Series in Statistics. Springer, Berlin, 2004.
[10] T. R. Hughes, M. J. Marton, A. R. Jones, et al. Functional discovery via a compendium of expression proiles. Cell, 102:109–126, 2000.
[11] M. M. Van Hulle. Edgeworth approximation of multivariate diﬀerential entropy.
Neural Computation, 17(9):1903–1910, 2005.
[12] S. Khan, S. Bandyopadhyay, A. Ganguly, and S. Saigal. Relative performance of
mutual information estimation methods for quantifying the dependence among shortand noisy data. Physical Review E, 76:026209, 2007.
[13] A. Kraskov, H. Stögbauer, and P. Grassberger. Estimating mutual information.
Physical Review E, 69:066138, 2004.
[14] M. S. Pepe. Evaluation of Medical Tests for Classification and Prediction. Oxford
Press, 2003.
[15] I. Priness, O. Maimon, and I. Ben-Gal. Evaluation of gene-expression clustering via
mutual information distance measure. BMC Bioinformatics, 8:111, 2007.
[16] Y. Saeys, I. Inza, and P. Larranaga. A review of feature selection techniques inbioinformatics. Bioinformatics, 23(19):2507–2517, 2007.
[17] B. T. Schlitt, K. Palin, J. Rung, S. Dietmann, M. Lappe, E. Ukkonen, and A. Brazma.
From gene networks to gene function. Genome Research, 13:2568–2576, 2003.
[18] B. Schölkopf and A. J. Smola. Learning with Kernels. MIT Press, Cambridge, MA,2002.

Mutual Information Estimation Reveals Global Associations

19[19] E. Segal, N. Friedman, D. Koller, and A. Regev. A module map showing conditional
activity of expression modules in cancer. Nature Genetics, 36(10):1090–1098, 2004.
[20] B. W. Silverman. Density Estimation for Statistics and Data Analysis. Chapman &Hall/CRC, April 1986.
[21] A. Subramanian, P. Tamayo, V. K. Mootha, S. Mukherjee, B. L. Ebert, M. A.
Gillette, A. Paulovich, S. L. Pomeroy, T. R. Golub, E. S. Lander, and J. P. Mesirov.
Gene set enrichment analysis: A knowledge-based approach for interpreting genomewide expression proﬁles. Proceedings of the National Academy of Sciences of theUSA, 102:15545–15550, 2005.
[22] K. Torkkola. Feature extraction by non-parametric mutual information maximization. Journal of Machine Learning Research, 3:1415–1438, 2003.
[23] S. van de Geer. Estimating a regression function. The Annals of Statistics, 18(2):907–
924, 1990.
[24] S. van de Geer. Empirical Processes in M-Estimation. Cambridge University Press,
2000.
[25] A. W. van der Vaart and J. A. Wellner. Weak Convergence and Empirical Processes.
With Applications to Statistics. Springer, New York, 1996.