# Multilevel Modeling with Structured Penalties for Classification from
Imaging Genetics data

Pascal Lu${}^{1,2(\text{\Letter})}$ 1\. Sorbonne Universités, UPMC Université
Paris 06, Inserm, CNRS, Institut du cerveau et la moelle (ICM), AP-HP -
Hôpital Pitié-Salpêtrière,
Boulevard de l’hôpital, 75013, Paris, France
<EMAIL_ADDRESS>
2\. INRIA Paris, ARAMIS project-team, 75013, Paris, France
<EMAIL_ADDRESS>Olivier Colliot1,2 1\. Sorbonne Universités, UPMC
Université Paris 06, Inserm, CNRS, Institut du cerveau et la moelle (ICM), AP-
HP - Hôpital Pitié-Salpêtrière,
Boulevard de l’hôpital, 75013, Paris, France
<EMAIL_ADDRESS>
2\. INRIA Paris, ARAMIS project-team, 75013, Paris, France
<EMAIL_ADDRESS>and
the Alzheimer’s Disease Neuroimaging Initiative  1\. Sorbonne Universités,
UPMC Université Paris 06, Inserm, CNRS, Institut du cerveau et la moelle
(ICM), AP-HP - Hôpital Pitié-Salpêtrière,
Boulevard de l’hôpital, 75013, Paris, France
<EMAIL_ADDRESS>
2\. INRIA Paris, ARAMIS project-team, 75013, Paris, France
<EMAIL_ADDRESS>
###### Abstract

In this paper, we propose a framework for automatic classification of patients
from multimodal genetic and brain imaging data by optimally combining them.
Additive models with unadapted penalties (such as the classical group lasso
penalty or $\displaystyle\ell_{1}$-multiple kernel learning) treat all
modalities in the same manner and can result in undesirable elimination of
specific modalities when their contributions are unbalanced. To overcome this
limitation, we introduce a multilevel model that combines imaging and genetics
and that considers joint effects between these two modalities for diagnosis
prediction. Furthermore, we propose a framework allowing to combine several
penalties taking into account the structure of the different types of data,
such as a group lasso penalty over the genetic modality and a
$\displaystyle\ell_{2}$-penalty on imaging modalities. Finally, we propose a
fast optimization algorithm, based on a proximal gradient method. The model
has been evaluated on genetic (single nucleotide polymorphisms - SNP) and
imaging (anatomical MRI measures) data from the ADNI database, and compared to
additive models [13, 15]. It exhibits good performances in AD diagnosis; and
at the same time, reveals relationships between genes, brain regions and the
disease status.

##  1 Introduction

The research area of imaging genetics studies the association between genetic
and brain imaging data [8]. A large number of papers studied the relationship
between genetic and neuroimaging data by considering that a phenotype can be
explained by a sum of effects from genetic variants. These multivariate
approaches use partial least squares [16], sparse canonical correlation
analysis [17], sparse regularized linear regression with a
$\displaystyle\ell_{1}$-penalty [10], group lasso penalty [12, 11], or
Bayesian model that links genetic variants to imaging regions and imaging
regions to the disease status [9].

But another interesting problem is about combining genetic and neuroimaging
data for automatic classification of patients. In particular, machine learning
methods have been used to build predictors for heterogeneous data, coming from
different modalities for brain disease diagnosis, such as Alzheimer’s disease
(AD) diagnosis. However, challenging issues are high-dimensional data, small
number of observations, the heterogeneous nature of data, and the weight for
each modality.

A framework that is commonly used to combine heterogeneous data is multiple
kernel learning (MKL) [6]. In MKL, each modality is represented by a kernel
(usually a linear kernel). The decision function and weights for the kernel
are simultaneously learnt. Moreover, the group lasso [2, 3] is a way to
integrate structure inside data. However, the standard
$\displaystyle\ell_{1}$-MKL and group lasso may eliminate modalities that have
a weak contribution. In particular, for AD, imaging data already provides good
results for its diagnosis. To overcome this problem, different papers have
proposed to use a $\displaystyle\ell_{1,p}$-penalty [7] to combine optimally
different modalities [13, 14].

These approaches do not consider potential effects between genetic and imaging
data for diagnosis prediction, as they only capture brain regions and SNPs
separately taken. Moreover, they put on the same level genetic and imaging
data, although these data do not provide the same type of information: given
only APOE genotyping, subjects can be classified according to their risk to
develop AD in the future; on the contrary, imaging data provides a photography
of the subject’s state at the present time.

Thereby, we propose a new framework that makes hierarchical the parameters and
considers interactions between genetic and imaging data for AD diagnosis. We
started with the idea that learning AD diagnosis from imaging data already
provides good results. Then, we considered that the decision function
parameters learnt from imaging data could be modulated, depending on each
subject’s genetic data. In other words, genes would express themselves through
these parameters. Considering a linear regression that links these parameters
and the genetic data, it leads to a multilevel model between imaging and
genetics. Our method also proposes potential relations between genetic and
imaging variables, if both of them are simultaneously related to AD. This
approach is different from the modeling proposed by [9], where imaging
variables are predicted from genetic variables, and diagnosis is predicted
from imaging variables.

Furthermore, current approaches [13, 14, 15] do not exploit data structure
inside each modality, as it is logical to group SNPs by genes, to expect
sparsity between genes (all genes are not linked to AD) and to enforce a
smooth regularization over brain regions for imaging modality. Thus, we have
imposed specific penalties for each modality by using a
$\displaystyle\ell_{2}$-penalty on the imaging modality, and a group lasso
penalty over the genetic modality. It models the mapping of variants into
genes, providing a better understanding of the role of genes in AD.

To learn all the decision function parameters, a fast optimization algorithm,
based on a proximal gradient method, has been developed. Finally, we have
evaluated our model on 1,107 genetic (SNP) and 114 imaging (anatomical MRI
measures) variables from the ADNI database111http://adni.loni.usc.edu and
compared it to additive models [13, 15].

##  2 Model set-up

###  2.1 Multilevel Logistic Regression with Structured Penalties

Let
$\displaystyle\\{(\mathbf{x}_{\mathcal{G}}^{k},\mathbf{x}_{\mathcal{I}}^{k},y^{k}),k=1,\ldots,N\\}$
be a set of labeled data, with
$\displaystyle\mathbf{x}_{\mathcal{G}}^{k}\in\mathbb{R}^{|\mathcal{G}|}$
(genetic data), and
$\displaystyle\mathbf{x}_{\mathcal{I}}^{k}\in\mathbb{R}^{|\mathcal{I}|}$
(imaging data) and $\displaystyle y^{k}\in\\{0,1\\}$ (diagnosis). Genetic,
imaging and genetic-imaging cross products training data are assumed centered
and normalized.

We propose the following Multilevel Logistic Regression model:

$\displaystyle
p(y=1|\mathbf{x}_{\mathcal{G}},\mathbf{x}_{\mathcal{I}})=\sigma\left(\boldsymbol{\alpha}(\mathbf{x}_{\mathcal{G}})^{\top}\mathbf{x}_{\mathcal{I}}+\alpha_{0}(\mathbf{x}_{\mathcal{G}})\right)\quad\text{with}~{}\sigma:x\mapsto\frac{1}{1+\text{e}^{-x}}$

where $\displaystyle\alpha_{0}(\mathbf{x}_{\mathcal{G}})$ is the intercept and
$\displaystyle\boldsymbol{\alpha}(\mathbf{x}_{\mathcal{G}})\in\mathbb{R}^{|\mathcal{I}|}$
is the parameter vector. On the contrary of the classical logistic regression
model, we propose a multilevel model, for which the parameter vector
$\displaystyle\boldsymbol{\alpha}(\mathbf{x}_{\mathcal{G}})$ and the intercept
$\displaystyle\alpha_{0}(\mathbf{x}_{\mathcal{G}})$ depend on genetic data
$\displaystyle\mathbf{x}_{\mathcal{G}}$.

This is to be compared to an additive model, where the diagnosis is directly
deduced from genetic and imaging data put at the same level. We assume that
$\displaystyle\boldsymbol{\alpha}$ and $\displaystyle\alpha_{0}$ are affine
functions of genetic data $\displaystyle\mathbf{x}_{\mathcal{G}}$:

$\displaystyle\boldsymbol{\alpha}(\mathbf{x}_{\mathcal{G}})=\mathbf{W}\mathbf{x}_{\mathcal{G}}+\boldsymbol{\beta}_{\mathcal{I}}\quad\text{and}\quad\alpha_{0}(\mathbf{x}_{\mathcal{G}})=\boldsymbol{\beta}_{\mathcal{G}}^{\top}\mathbf{x}_{\mathcal{G}}+\beta_{0}$

where
$\displaystyle\mathbf{W}\in\mathcal{M}_{|\mathcal{I}|,|\mathcal{G}|}(\mathbb{R})$,
$\displaystyle\boldsymbol{\beta}_{\mathcal{I}}\in\mathbb{R}^{|\mathcal{I}|}$,
$\displaystyle\boldsymbol{\beta}_{\mathcal{G}}\in\mathbb{R}^{|\mathcal{G}|}$
and $\displaystyle\beta_{0}\in\mathbb{R}$. Therefore, the probability becomes
$\displaystyle
p(y=1|\mathbf{x}_{\mathcal{G}},\mathbf{x}_{\mathcal{I}})=\sigma\left(\mathbf{x}_{\mathcal{G}}^{\top}\mathbf{W}^{\top}\mathbf{x}_{\mathcal{I}}+\boldsymbol{\beta}_{\mathcal{I}}^{\top}\mathbf{x}_{\mathcal{I}}+\boldsymbol{\beta}_{\mathcal{G}}^{\top}\mathbf{x}_{\mathcal{G}}+\beta_{0}\right)$.
Figure 1 summarizes the relations between parameters.

$\displaystyle
y$$\displaystyle\boldsymbol{\alpha}(\mathbf{x}_{\mathcal{G}})$$\displaystyle\mathbf{x}_{\mathcal{I}}$$\displaystyle\alpha_{0}(\mathbf{x}_{\mathcal{G}})$$\displaystyle\beta_{0}$$\displaystyle\boldsymbol{\beta}_{\mathcal{G}}$$\displaystyle\mathbf{x}_{\mathcal{G}}$$\displaystyle\boldsymbol{\beta}_{\mathcal{I}}$$\displaystyle\mathbf{W}$

Figure 1: The disease status $\displaystyle y$ is predicted from imaging data
$\displaystyle\mathbf{x}_{\mathcal{I}}$ and the parameters
$\displaystyle\beta_{0}(\mathbf{x}_{\mathcal{G}}),\boldsymbol{\beta}(\mathbf{x}_{\mathcal{G}})$
(which are computed from genetic data $\displaystyle\mathbf{x}_{\mathcal{G}}$)

The parameters
$\displaystyle\mathbf{W},\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}},\beta_{0}$
are obtained by minimizing the objective:

$\displaystyle\displaystyle
S(\mathbf{W},\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}},\beta_{0})$
$\displaystyle\displaystyle=R_{N}(\mathbf{W},\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}},\beta_{0})+\Omega(\mathbf{W},\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}})$
$\displaystyle\displaystyle\text{with
}R_{N}(\mathbf{W},\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}},\beta_{0})$
$\displaystyle\displaystyle=\frac{1}{N}\sum_{k=1}^{N}\left\\{-y^{k}\left({(\mathbf{x}_{\mathcal{G}}^{k})}^{\top}\mathbf{W}^{\top}\mathbf{x}_{\mathcal{I}}^{k}+\boldsymbol{\beta}_{\mathcal{I}}^{\top}\mathbf{x}_{\mathcal{I}}^{k}+\boldsymbol{\beta}_{\mathcal{G}}^{\top}\mathbf{x}_{\mathcal{G}}^{k}+\beta_{0}\right)\right.$
$\displaystyle\displaystyle\quad\quad\quad\quad\quad\quad\left.+\log\left(1+\text{e}^{(\mathbf{x}_{\mathcal{G}}^{k})^{\top}\mathbf{W}^{\top}\mathbf{x}_{\mathcal{I}}^{k}+\boldsymbol{\beta}_{\mathcal{I}}^{\top}\mathbf{x}_{\mathcal{I}}^{k}+\boldsymbol{\beta}_{\mathcal{G}}^{\top}\mathbf{x}_{\mathcal{G}}^{k}+\beta_{0}}\right)\right\\}$
$\displaystyle\displaystyle\text{and
}\quad\quad\Omega(\mathbf{W},\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}})$
$\displaystyle\displaystyle=\lambda_{W}\Omega_{W}(\mathbf{W})+\lambda_{\mathcal{I}}\Omega_{\mathcal{I}}(\boldsymbol{\beta}_{\mathcal{I}})+\lambda_{\mathcal{G}}\Omega_{\mathcal{G}}(\boldsymbol{\beta}_{\mathcal{G}})$

$\displaystyle\Omega_{W}$, $\displaystyle\Omega_{\mathcal{I}}$,
$\displaystyle\Omega_{\mathcal{G}}$ are respectively the penalties for
$\displaystyle\mathbf{W}$, $\displaystyle\boldsymbol{\beta}_{\mathcal{I}}$,
$\displaystyle\boldsymbol{\beta}_{\mathcal{G}}$, whereas
$\displaystyle\lambda_{W}>0$, $\displaystyle\lambda_{\mathcal{I}}>0$,
$\displaystyle\lambda_{\mathcal{G}}>0$ are respectively the regularization
parameters for $\displaystyle\Omega_{W}$, $\displaystyle\Omega_{\mathcal{I}}$,
$\displaystyle\Omega_{\mathcal{G}}$.

Genetic data are a sequence of single-polymorphism nucleotides (SNP) counted
by minor allele. A SNP can belong (or not) to one gene $\displaystyle{\ell}$
(or more) and therefore participate in the production of proteins that
interact inside pathways. We decided to group SNPs by genes, and designed a
penalty to enforce sparsity between genes and regularity inside genes. Given
that some SNPs may belong to multiple genes, the group lasso with overlap
penalty [4] is more suitable, with genes as groups. To deal with this penalty,
an overlap expansion is performed. Given
$\displaystyle\mathbf{x}\in\mathbb{R}^{|\mathcal{G}|}$ a subject’s feature
vector, a new feature vector is created
$\displaystyle\widetilde{\mathbf{x}}=\left(\mathbf{x}_{\mathcal{G}_{1}}^{\top},\ldots,\mathbf{x}_{\mathcal{G}_{L}}^{\top}\right)^{\top}\in\mathbb{R}^{\sum_{{\ell}=1}^{L}|\mathcal{G}_{\ell}|}$,
defined by the concatenation of copies of the genetic data restricted by group
$\displaystyle\mathcal{G}_{\ell}$. Similarly, the same expansion is performed
on $\displaystyle{\boldsymbol{\beta}}_{\mathcal{G}},\mathbf{W}$ to obtain
$\displaystyle\widetilde{\boldsymbol{\beta}}_{\mathcal{G}}\in\mathbb{R}^{\sum_{{\ell}=1}^{L}|\mathcal{G}_{\ell}|}$
and
$\displaystyle\widetilde{\mathbf{W}}\in\mathbb{R}^{|\mathcal{I}|\times(\sum_{{\ell}=1}^{L}|\mathcal{G}_{\ell}|)}$.
This group lasso with overlap penalty is used for the matrix
$\displaystyle\mathbf{W}$ and for
$\displaystyle{\boldsymbol{\beta}}_{\mathcal{G}}$.

For imaging variables, the ridge penalty is considered:
$\displaystyle\Omega_{\mathcal{I}}(\boldsymbol{\beta}_{\mathcal{I}})=\left\|\boldsymbol{\beta}_{\mathcal{I}}\right\|_{2}^{2}$.
In particular, brain diseases usually have a diffuse anatomical pattern of
alteration throughout the brain and therefore, regularity is usually required
for the imaging parameter. Finally, $\displaystyle\Omega$ is defined by:

$\displaystyle\Omega\left(\widetilde{\mathbf{W}},\widetilde{\boldsymbol{\beta}}_{\mathcal{G}},\boldsymbol{\beta}_{\mathcal{I}}\right)=\lambda_{W}\sum_{i=1}^{|\mathcal{I}|}\sum_{\ell=1}^{L}\theta_{\mathcal{G}_{\ell}}\left\|\widetilde{\mathbf{W}}_{i,\mathcal{G}_{\ell}}\right\|_{2}+\lambda_{\mathcal{I}}\left\|\widetilde{\boldsymbol{\beta}}_{\mathcal{I}}\right\|_{2}+\lambda_{\mathcal{G}}\sum_{{\ell}=1}^{L}\theta_{\mathcal{G}_{\ell}}\left\|\widetilde{\boldsymbol{\beta}}_{\mathcal{G}_{\ell}}\right\|_{2}$

###  2.2 Minimization of $\displaystyle
S(\mathbf{W},\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{I}},\beta_{0})$

From now on, and for simplicity reasons,
$\displaystyle\widetilde{\mathbf{W}}$,
$\displaystyle\widetilde{\boldsymbol{\beta}}$ and
$\displaystyle\widetilde{\mathbf{x}}$ are respectively denoted as
$\displaystyle\mathbf{W}$, $\displaystyle\boldsymbol{\beta}$ and
$\displaystyle{\mathbf{x}}$. Let $\displaystyle\Phi$ be the function that
reshapes a matrix of
$\displaystyle\mathcal{M}_{|\mathcal{I}|,|\mathcal{G}|}(\mathbb{R})$ to a
vector of $\displaystyle\mathbb{R}^{|\mathcal{I}|\times|\mathcal{G}|}$ (i.e.
$\displaystyle\mathbf{W}_{i,g}=\Phi(\mathbf{W})_{i|\mathcal{G}|+g}$):

$\displaystyle\Phi:\mathbf{W}\mapsto((\mathbf{W}_{1,1},\ldots,\mathbf{W}_{1,|\mathcal{G}|}),\ldots,(\mathbf{W}_{|\mathcal{I}|,1},\ldots,\mathbf{W}_{|\mathcal{I}|,|\mathcal{G}|}))$

We will estimate $\displaystyle\Phi(\mathbf{W})$ and then reshape it to obtain
$\displaystyle\mathbf{W}$. The algorithm developed is based on a proximal
gradient method [1, 5].

The parameters
$\displaystyle\mathbf{w}^{(t+1)}=\left(\Phi\left(\mathbf{W}^{(t+1)}\right),\boldsymbol{\beta}_{\mathcal{I}}^{(t+1)},\boldsymbol{\beta}_{\mathcal{G}}^{(t+1)},{\beta}_{0}^{(t+1)}\right)$
are updated with:

$\displaystyle\displaystyle\mathbf{w}^{(t+1)}$
$\displaystyle\displaystyle=\mathop{\text{argmin}}_{\mathbf{w}}R_{N}\left(\mathbf{w}\right)+\left[\mathbf{w}-\mathbf{w}^{(t)}\right]^{\top}\nabla
R_{N}\left(\mathbf{w}^{(t)}\right)+\frac{1}{2\varepsilon}\left\|\mathbf{w}-\mathbf{w}^{(t)}\right\|_{2}^{2}+\Omega(\mathbf{w})$
$\displaystyle\displaystyle=\mathop{\text{argmin}}_{\mathbf{w}}\left\\{\frac{1}{2}\left\|\boldsymbol{\omega}^{(t)}-\mathbf{w}^{(t)}\right\|_{2}^{2}+\varepsilon\Omega(\mathbf{w})\right\\}\text{
with }\boldsymbol{\omega}^{(t)}=\mathbf{w}^{(t)}-\varepsilon\nabla
R_{N}\left(\mathbf{w}^{(t)}\right)$

The idea is to update $\displaystyle\mathbf{w}^{(t+1)}$ from
$\displaystyle\mathbf{w}^{(t)}$ with a Newton-type algorithm without the
constraint $\displaystyle\Omega$ given a stepsize $\displaystyle\varepsilon$,
and then to project the result onto the compact set defined by
$\displaystyle\Omega$. Regarding the stepsize $\displaystyle\varepsilon$, a
backtracking line search [5] is performed. Let
$\displaystyle\widehat{G}\left(\mathbf{w}^{(t)},\varepsilon\right)=\frac{1}{\varepsilon}\left[\mathbf{w}^{(t)}-\mathbf{w}^{(t+1)}\right]$
be the step in the proximal gradient update. A line search is performed over
$\displaystyle\varepsilon$ until the inequality is reached:

$\displaystyle\displaystyle R_{N}\left(\mathbf{w}^{(t+1)}\right)$
$\displaystyle\displaystyle\leq
R_{N}\left(\mathbf{w}^{(t)}\right)-\varepsilon\nabla
R_{N}\left(\mathbf{w}^{(t)}\right)^{\top}\widehat{G}\left(\mathbf{w}^{(t)},\varepsilon\right)+\frac{\varepsilon}{2}\left\|\widehat{G}\left(\mathbf{w}^{(t)},\varepsilon\right)\right\|_{2}^{2}$

The minimization algorithm stops when
$\displaystyle\left|S\left(\mathbf{w}^{(t+1)}\right)-S\left(\mathbf{w}^{(t)}\right)\right|\leq\eta\left|S\left(\mathbf{w}^{(t)}\right)\right|$,
where $\displaystyle\eta=10^{-5}$. The whole algorithm is summarized below:

1

2Input:
$\displaystyle\\{(\mathbf{x}_{\mathcal{I}}^{k},\mathbf{x}_{\mathcal{G}}^{k},y^{k}),k=1,\ldots,N\\}$,
$\displaystyle\delta=0.8$, $\displaystyle\varepsilon_{0}=1$,
$\displaystyle\eta=10^{-5}$ ;

3 Initialization: $\displaystyle\mathbf{W}=\mathbf{0}$,
$\displaystyle\boldsymbol{\beta}_{\mathcal{I}}=\mathbf{0}$,
$\displaystyle\boldsymbol{\beta}_{\mathcal{G}}=\mathbf{0}$,
$\displaystyle\beta_{0}=0$ and continue = True ;

4 while  _continue_ do

5 $\displaystyle\varepsilon=\varepsilon_{0}$;

6

7 $\displaystyle
R_{N}=R_{N}\left({\mathbf{W}},{\boldsymbol{\beta}}_{\mathcal{I}},{\boldsymbol{\beta}}_{\mathcal{G}},{\beta}_{0}\right)$;

8

9 $\displaystyle\displaystyle{\nabla
R_{N}=\frac{1}{N}\sum_{k=1}^{N}\begin{pmatrix}\Phi\left({(\mathbf{x}_{\mathcal{I}}^{k})}^{\top}\mathbf{x}_{\mathcal{G}}^{k}\right)\\\
\mathbf{x}_{\mathcal{I}}^{k}\\\ \mathbf{x}_{\mathcal{G}}^{k}\\\
1\end{pmatrix}\left[\sigma\left({(\mathbf{x}_{\mathcal{G}}^{k})}^{\top}\mathbf{W}^{\top}\mathbf{x}_{\mathcal{I}}^{k}+\boldsymbol{\beta}_{\mathcal{I}}^{\top}\mathbf{x}_{\mathcal{I}}^{k}+\boldsymbol{\beta}_{\mathcal{G}}^{\top}\mathbf{x}_{\mathcal{G}}^{k}+\beta_{0}\right)-y^{k}\right]}$

10
$\displaystyle\left(\widehat{\mathbf{W}},\widehat{\boldsymbol{\beta}}_{\mathcal{I}},\widehat{\boldsymbol{\beta}}_{\mathcal{G}},\widehat{\beta}_{0},\widehat{G}\right)=$
Algo_2$\displaystyle({\mathbf{W}},{\boldsymbol{\beta}}_{\mathcal{I}},{\boldsymbol{\beta}}_{\mathcal{G}},{\beta}_{0},\nabla
R_{N},\varepsilon)$;

11

12 while  _$\displaystyle
R_{N}\left(\widehat{\mathbf{W}},\widehat{\boldsymbol{\beta}}_{\mathcal{I}},\widehat{\boldsymbol{\beta}}_{\mathcal{G}},\widehat{\beta}_{0}\right)
>R_{N}-\varepsilon\nabla
R_{N}^{\top}\widehat{G}+\frac{\varepsilon}{2}\|\widehat{G}\|_{2}^{2}$_ do

13 $\displaystyle\varepsilon=\delta\varepsilon$ and
$\displaystyle\left(\widehat{\mathbf{W}},\widehat{\boldsymbol{\beta}}_{\mathcal{I}},\widehat{\boldsymbol{\beta}}_{\mathcal{G}},\widehat{\beta}_{0},\widehat{G}\right)=$
Algo_2$\displaystyle({\mathbf{W}},{\boldsymbol{\beta}}_{\mathcal{I}},{\boldsymbol{\beta}}_{\mathcal{G}},{\beta}_{0},\nabla
R_{N},\varepsilon)$;

14

15 end while

16 continue =
$\displaystyle\left|S\left(\widehat{\mathbf{W}},\widehat{\boldsymbol{\beta}}_{\mathcal{I}},\widehat{\boldsymbol{\beta}}_{\mathcal{G}},\widehat{\beta}_{0}\right)-S\left(\mathbf{W},{\boldsymbol{\beta}}_{\mathcal{I}},{\boldsymbol{\beta}}_{\mathcal{G}},\beta_{0}\right)\right|\mathop{>}^{\text{?}}\eta\left|S\left(\mathbf{W},{\boldsymbol{\beta}}_{\mathcal{I}},{\boldsymbol{\beta}}_{\mathcal{G}},\beta_{0}\right)\right|$

17 $\displaystyle{\mathbf{W}}=\widehat{\mathbf{W}}$,
$\displaystyle{\boldsymbol{\beta}}_{\mathcal{I}}=\widehat{\boldsymbol{\beta}}_{\mathcal{I}}$,
$\displaystyle{\boldsymbol{\beta}}_{\mathcal{G}}=\widehat{\boldsymbol{\beta}}_{\mathcal{G}}$,
$\displaystyle\beta_{0}=\widehat{\beta}_{0}$;

18

19 end while

return
$\displaystyle\left(\mathbf{W},{\boldsymbol{\beta}}_{\mathcal{I}},{\boldsymbol{\beta}}_{\mathcal{G}},\beta_{0}\right)$

Algorithm 1 Training the multilevel logistic regression

1 Input:
$\displaystyle({\mathbf{W}},{\boldsymbol{\beta}}_{\mathcal{I}},{\boldsymbol{\beta}}_{\mathcal{G}},{\beta}_{0})$
(parameters), $\displaystyle\nabla R_{N}$ (gradient),
$\displaystyle\varepsilon$ (stepsize) ;

2

3Compute
$\displaystyle\boldsymbol{\omega}=\boldsymbol{\beta}-\varepsilon\nabla_{({\mathbf{W}},{\boldsymbol{\beta}}_{\mathcal{I}},{\boldsymbol{\beta}}_{\mathcal{G}})}R_{N}$
;

4

5Update
$\displaystyle\widehat{\mathbf{W}}_{\mathcal{G}_{\ell},i}=\max\left(0,1-\frac{\varepsilon\lambda_{\mathcal{G}}\theta_{\mathcal{G}_{\ell}}}{\left\|\boldsymbol{\omega}^{(t)}_{\mathcal{G}_{\ell}+i|\mathcal{G}|}\right\|_{2}}\right)\boldsymbol{\omega}_{\mathcal{G}_{\ell}+i|\mathcal{G}|}$
for $\displaystyle(i,\ell)\in\llbracket
1,|\mathcal{I}|\rrbracket\times\llbracket 1,L\rrbracket$ ;

6

7Update
$\displaystyle\widehat{\boldsymbol{\beta}}_{\mathcal{I}}=\frac{\boldsymbol{\omega}_{\mathcal{I}+|\mathcal{G}||\mathcal{I}|}}{1+2\varepsilon\lambda_{\mathcal{I}}}$
(imaging modality) ;

8

9Update
$\displaystyle\widehat{\boldsymbol{\beta}}_{\mathcal{G}_{\ell}}=\max\left(0,1-\frac{\varepsilon\lambda_{\mathcal{G}}\theta_{\mathcal{G}_{\ell}}}{\left\|\boldsymbol{\omega}^{(t)}_{\mathcal{G}_{\ell}+(|\mathcal{G}|+1)|\mathcal{I}|}\right\|_{2}}\right)\boldsymbol{\omega}_{\mathcal{G}_{\ell}+(|\mathcal{G}|+1)|\mathcal{I}|}$
for $\displaystyle\ell\in\llbracket 1,L\rrbracket$ ;

10

11Update $\displaystyle\widehat{\beta}_{0}=\beta_{0}-\varepsilon\frac{\partial
R_{N}}{\partial\beta_{0}}$ and
$\displaystyle\widehat{G}=\frac{1}{\varepsilon}\left[\begin{pmatrix}\Phi({\mathbf{W}})\\\
{\boldsymbol{\beta}}_{\mathcal{I}}\\\ {\boldsymbol{\beta}}_{\mathcal{G}}\\\
{{\beta}}_{0}\end{pmatrix}-\begin{pmatrix}\Phi(\widehat{\mathbf{W}})\\\
\widehat{\boldsymbol{\beta}}_{\mathcal{I}}\\\
\widehat{\boldsymbol{\beta}}_{\mathcal{G}}\\\
\widehat{{\beta}}_{0}\end{pmatrix}\right]$ ;

return
$\displaystyle\left(\widehat{\mathbf{W}},\widehat{\boldsymbol{\beta}}_{\mathcal{I}},\widehat{\boldsymbol{\beta}}_{\mathcal{G}},\widehat{\beta}_{0},\widehat{G}\right)$

Algorithm 2 Parameter update

##  3 Experimental results

###  3.1 Dataset

The ADNI1 GWAS dataset from ADNI studied 707 subjects, with 156 Alzheimer’s
Disease patients (denoted AD), 196 MCI patients at baseline who progressed to
AD (denoted pMCI, as progressive MCI), 150 MCI patients who remain stable
(denoted sMCI, as stable MCI) and 201 healthy control subjects (denoted CN).

In ADNI1 GWAS dataset, 620,901 SNPs have been genotyped, but we selected 1,107
SNPs based on the 44 first top genes related to AD (from
AlzGene222http://www.alzgene.org) and on the Illumina annotation using the
Genome build 36.2. Group weighting for genes is based on gene size: for group
$\displaystyle\mathcal{G}_{\ell}$, the weight
$\displaystyle\theta_{\mathcal{G}_{\ell}}=\sqrt{|\mathcal{G}_{\ell}|}$ ensures
that the penalty term is of the order of the number of parameters of the
group.

The parameter $\displaystyle\lambda_{\mathcal{G}}$ influences the number of
groups that are selected by the model. In particular, the group
$\displaystyle\mathcal{G}_{\ell}$ enters in the model during the first
iteration if
$\displaystyle\left\|\nabla_{\boldsymbol{\beta}_{\mathcal{G}_{\ell}}}R_{N}(\mathbf{0})\right\|_{2}>\lambda_{\mathcal{G}}\theta_{\mathcal{G}_{\ell}}$.
This inequality gives an upper bound for $\displaystyle\lambda_{\mathcal{G}}$.
The same remark can be done for $\displaystyle\lambda_{W}$. Regarding MRI
modality, we used the segmentation of FreeSurfer which gives the volume of
subcortical regions (44 features) and the average cortical region thickness
(70 features). Therefore, there are $\displaystyle 1,107\times 114=126,198$
parameters to infer for $\displaystyle\mathbf{W}$, $\displaystyle 114$
parameters for $\displaystyle\boldsymbol{\beta}_{\mathcal{I}}$ and
$\displaystyle 1,107$ parameters for
$\displaystyle\boldsymbol{\beta}_{\mathcal{G}}$.

###  3.2 Results

We ran our multilevel model and compared it to the logistic regression applied
to one single modality with simple penalties (lasso, group lasso, ridge), to
additive models ([13], [15] EasyMKL with a linear kernel for each modality,
and the model $\displaystyle
p(y=1|\mathbf{x}_{\mathcal{G}},\mathbf{x}_{\mathcal{I}})=\sigma\left(\boldsymbol{\beta}_{\mathcal{I}}^{\top}\mathbf{x}_{\mathcal{I}}+\boldsymbol{\beta}_{\mathcal{G}}^{\top}\mathbf{x}_{\mathcal{G}}+\beta_{0}\right)$
with our algorithm under the constraint
$\displaystyle\boldsymbol{\beta}_{\mathcal{G}}\neq\mathbf{0}$), and to the
multiplicative model with $\displaystyle\mathbf{W}$ only, where $\displaystyle
p(y=1|\mathbf{x}_{\mathcal{G}},\mathbf{x}_{\mathcal{I}})=\sigma\left(\mathbf{x}_{\mathcal{G}}^{\top}\mathbf{W}^{\top}\mathbf{x}_{\mathcal{I}}+\beta_{0}\right)$.
We considered two classification tasks: “AD versus CN” and “pMCI versus CN”.
Four measures are used: the sensitivity (Sen), the specificity (Spe), the
precision (Pre) and the balanced accuracy between the sensitivity and the
specificity (BAcc). A 10-fold cross validation is performed. The parameters
$\displaystyle\lambda_{W},\lambda_{\mathcal{I}},\lambda_{\mathcal{G}}$ are
optimised between $\displaystyle[10^{-3},1]$. Classification results for these
tasks are shown on table 1. It typically takes between 5 and 8 minutes to
learn the parameters.

Table 1: Classification results for different modalities and methods

|  AD versus CN (%)
---|---
Modality |  Method & Penalty |  Sen |  Spe |  Pre | BAcc
SNPs only | logistic regression (lasso $\displaystyle\ell_{1}$) | 69.4 | 77.5 | 71.1 | 73.4
SNPs grouped by genes | logistic regression (group lasso) | 69.4 | 77.5 | 71.1 | 73.4
MRI (cortical) | logistic regression (ridge $\displaystyle\ell_{2}$) | 84.4 | 89.5 | 87.1 | 86.9
MRI (subcortical) | logistic regression (ridge $\displaystyle\ell_{2}$) | 80.0 | 86.0 | 83.2 | 83.0
SNP + MRI (all) |  [15] EasyMKL, _Aiolli et al._ | 89.4 | 85.0 | 83.0 | 87.2
SNP + MRI (all) |  [13] _Wang et al._ | 89.4 | 88.0 | 85.7 | 88.7
SNP + MRI (all) | additive model ($\displaystyle\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}}$ only) | 88.8 | 89.5 | 87.6 | 89.1
SNP + MRI (all) | multiplicative model ($\displaystyle\mathbf{W}$ only) | 89.4 | 87.0 | 85.0 | 88.2
SNP + MRI (all) | multilevel model (all) | 90.6 | 87.0 | 85.5 | 88.8

|  pMCI versus CN (%)
---|---
Modality |  Method & Penalty |  Sen |  Spe |  Pre | BAcc
SNPs only | logistic regression (lasso $\displaystyle\ell_{1}$) | 72.0 | 77.0 | 75.9 | 74.5
SNPs grouped by genes | logistic regression (group lasso) | 72.0 | 77.0 | 75.9 | 74.5
MRI (cortical) | logistic regression (ridge $\displaystyle\ell_{2}$) | 74.0 | 76.0 | 76.4 | 75.0
MRI (subcortical) | logistic regression (ridge $\displaystyle\ell_{2}$) | 73.0 | 76.5 | 76.6 | 74.7
SNP + MRI (all) |  [15] EasyMKL, _Aiolli et al._ | 77.0 | 73.5 | 75.1 | 75.3
SNP + MRI (all) |  [13] _Wang et al._ | 79.5 | 81.5 | 82.4 | 80.5
SNP + MRI (all) | additive model ($\displaystyle\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}}$ only) | 80.5 | 81.0 | 82.0 | 80.8
SNP + MRI (all) | multiplicative model ($\displaystyle\mathbf{W}$ only) | 81.0 | 81.5 | 82.9 | 81.3
SNP + MRI (all) | multilevel model (all) | 82.5 | 83.0 | 84.1 | 82.8

Regarding MRI features, the most important features (in weight) are the
left/right hippocampus, the left/right Amygdala, the left/right entorhinal and
the left middle temporal cortices. Regarding genetic features, the most
important features in weight are SNPs that belong to gene APOE (rs429358) for
both tasks “AD versus CN” and “pMCI versus CN”.

Figure 2: Overview of the reduced parameters
$\displaystyle\overline{\mathbf{W}}\in\mathbb{R}^{|\mathcal{I}|\times L}$,
$\displaystyle\overline{\boldsymbol{\beta}}_{\mathcal{I}}\in\mathbb{R}^{|\mathcal{I}|}$
and
$\displaystyle\overline{\boldsymbol{\beta}}_{\mathcal{G}}\in\mathbb{R}^{L}$
(learnt through the task “pMCI vs CN” for the whole model). For brain region
$\displaystyle i$ and gene $\displaystyle{\ell}$,
$\displaystyle\overline{\mathbf{W}}[i,{\ell}]=\max_{g\in\mathcal{G}_{\ell}}|\mathbf{W}[i,g]|$,
$\displaystyle\overline{\boldsymbol{\beta}}_{\mathcal{I}}[i]=|{\boldsymbol{\beta}}_{\mathcal{I}}[i]|$
and
$\displaystyle\overline{\boldsymbol{\beta}}_{\mathcal{G}}[\ell]=\max_{g\in\mathcal{G}_{\ell}}|{\boldsymbol{\beta}}_{\mathcal{G}}[g]|$.
Only some brain regions are shown in this figure.

Regarding the matrix $\displaystyle\mathbf{W}$, the couples (brain region,
gene) learnt through the task “pMCI versus CN” are shown on Fig. 2. It can be
seen that $\displaystyle\mathbf{W}$ has a sparse structure. Among the couples
(brain region, gene) that have non null coefficients for the both tasks “AD
versus CN” and “pMCI versus CN”, there are (Left Hippocampus, MGMT), (Right
Entorhinal, APOE) or (Left Middle Temporal, APOE). Only couples related to AD
are selected by the model.

We noticed that genes and brain regions strongly related to AD are captured by
the vectors $\displaystyle{\boldsymbol{\beta}}_{\mathcal{G}}$ and
$\displaystyle{\boldsymbol{\beta}}_{\mathcal{I}}$, whereas genes less strongly
related to AD are captured by the matrix $\displaystyle\mathbf{W}$. Coming
back to original formulation described in section 2.1, the contribution of the
function
$\displaystyle\alpha_{0}:\mathbf{x}_{\mathcal{G}}\mapsto{\boldsymbol{\beta}}_{\mathcal{G}}^{\top}\mathbf{x}_{\mathcal{G}}+\beta_{0}$
is much smaller (in terms of weights) than the function
$\displaystyle\boldsymbol{\alpha}:\mathbf{x}_{\mathcal{G}}\mapsto\mathbf{W}\mathbf{x}_{\mathcal{G}}+{\boldsymbol{\beta}}_{\mathcal{I}}$.
Furthermore, Fig. 2 shows that genetic data
$\displaystyle\mathbf{x}_{\mathcal{G}}$ tend to express through
$\displaystyle\mathbf{W}$, and thereby participate in the modulation of the
vector $\displaystyle{\boldsymbol{\alpha}}(\mathbf{x}_{\mathcal{G}})$.

We compared our approach to [13, 15], for which the codes are available. The
features that are selected by [13, 15] are similar to ours for each modality
taken separately. For instance, for [13] and the task “AD versus CN”, SNPs
that have the most important weights are in genes APOE (rs429358), BZW1
(rs3815501) and MGMT (rs7071424). However, the genetic parameter vector learnt
from [13] or [15] is not sparse, in contrary of ours. Furthermore, for [15],
the weight for the imaging kernel is nine times much larger than the weight
for the genetic kernel. These experiments show that the additive model with
adapted penalties for each modality provides better performances than [15],
but our additive, multiplicative and multilevel models provide similar
performances.

##  4 Conclusion

In this paper, we developed a novel approach to integrate genetic and brain
imaging data for prediction of disease status. Our multilevel model takes into
account potential interactions between genes and brain regions, but also the
structure of the different types of data though the use of specific penalties
within each modality. When applied to genetic and MRI data from the ADNI
database, the model was able to highlight brain regions and genes that have
been previously associated with AD, thereby demonstrating the potential of our
approach for imaging genetics studies in brain diseases.

Acknowledgments. The research leading to these results has received funding
from the program _Investissements d’avenir ANR-10-IAIHU-06_. We wish to thank
Theodoros Evgeniou for many useful insights.

##  Appendix 0.A Probabilistic formulation

This section proposes a probabilistic formulation for the model. The
conditional probability is given by $\displaystyle
p(y=1|\mathbf{x}_{\mathcal{G}},\mathbf{x}_{\mathcal{I}})=\sigma\left(\mathbf{x}_{\mathcal{G}}^{\top}\mathbf{W}^{\top}\mathbf{x}_{\mathcal{I}}+\boldsymbol{\beta}_{\mathcal{I}}^{\top}\mathbf{x}_{\mathcal{I}}+\boldsymbol{\beta}_{\mathcal{G}}^{\top}\mathbf{x}_{\mathcal{G}}+\beta_{0}\right)$.

  * –

For each region $\displaystyle i\in\mathcal{I}$ and gene
$\displaystyle\mathcal{G}_{\ell}$,
$\displaystyle\mathbf{W}_{i,\mathcal{G}_{\ell}}\sim\text{M-Laplace}(0,\lambda_{W})$
(M-Laplace stands for “Multi-Laplacian prior”). In other words:

$\displaystyle\displaystyle
p(\mathbf{W};\lambda_{W},\mathcal{G},\boldsymbol{\theta}_{\mathcal{G}})$
$\displaystyle\displaystyle\propto\prod_{i=1}^{|\mathcal{I}|}\prod_{\ell=1}^{L}\text{e}^{-\lambda_{W}{\theta}_{\mathcal{G}_{\ell}}\|\mathbf{W}_{i,\mathcal{G}_{\ell}}\|_{2}}$

  * –

For each region $\displaystyle i\in\mathcal{I}$,
$\displaystyle\boldsymbol{\beta}_{i}\sim\mathcal{N}\left(0,\frac{1}{2\lambda_{\mathcal{I}}}\right)$,
i.e. $\displaystyle
p(\boldsymbol{\beta}_{\mathcal{I}};\lambda_{\mathcal{I}})\propto\text{e}^{-\lambda_{\mathcal{I}}\|\boldsymbol{\beta}_{\mathcal{I}}\|_{2}^{2}}$

  * –

For each gene $\displaystyle\mathcal{G}_{\ell}$,
$\displaystyle\boldsymbol{\beta}_{\mathcal{G}_{\ell}}\sim\text{M-Laplace}(0,\lambda_{\mathcal{G}})$,
i.e.

$\displaystyle
p(\boldsymbol{\beta}_{\mathcal{G}};\lambda_{\mathcal{G}},\mathcal{G},\boldsymbol{\theta}_{\mathcal{G}})\propto\prod_{\ell=1}^{L}\text{e}^{-\lambda_{\mathcal{G}}{\theta}_{\mathcal{G}_{\ell}}\|{\boldsymbol{\beta}}_{\mathcal{G}_{\ell}}\|_{2}}$

Let $\displaystyle Y=(y^{1},\ldots,y^{N})$, $\displaystyle
X_{\mathcal{I}}=(\mathbf{x}_{\mathcal{I}}^{1},\ldots,\mathbf{x}_{\mathcal{I}}^{N})$
and $\displaystyle
X_{\mathcal{G}}=(\mathbf{x}_{\mathcal{G}}^{1},\ldots,\mathbf{x}_{\mathcal{G}}^{N})$.

The generative model is given by:

$\displaystyle
p(\mathbf{W},\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}},\beta_{0},Y,X_{\mathcal{I}},X_{\mathcal{G}};\lambda_{W},\lambda_{\mathcal{I}},\lambda_{\mathcal{G}},\mathcal{G},\boldsymbol{\theta}_{\mathcal{G}})$

$\displaystyle\displaystyle\mathop{=}^{\text{Bayes}}p(Y,X_{\mathcal{I}},X_{\mathcal{G}}|\mathbf{W},\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}})p(\mathbf{W};\lambda_{W},\mathcal{G},\boldsymbol{\theta}_{\mathcal{G}})p(\boldsymbol{\beta}_{\mathcal{I}};\lambda_{\mathcal{I}})p(\boldsymbol{\beta}_{\mathcal{G}};\lambda_{\mathcal{G}},\mathcal{G},\boldsymbol{\theta}_{\mathcal{G}})p(\beta_{0})$
$\displaystyle\displaystyle\mathop{=}^{\text{obs
iid}}\left(\prod_{k=1}^{N}p(y=y^{k},\mathbf{x}_{\mathcal{I}}^{k},\mathbf{x}_{\mathcal{G}}^{k}|\mathbf{W},\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}})\right)$
$\displaystyle\displaystyle\hskip
28.45274ptp(\mathbf{W};\lambda_{W},\mathcal{G},\boldsymbol{\theta}_{\mathcal{G}})p(\boldsymbol{\beta}_{\mathcal{I}};\lambda_{\mathcal{I}})p(\boldsymbol{\beta}_{\mathcal{G}};\lambda_{\mathcal{G}},\mathcal{G},\boldsymbol{\theta}_{\mathcal{G}})p(\beta_{0})$
$\displaystyle\displaystyle\propto\prod_{k=1}^{N}\sigma\left({(\mathbf{x}_{\mathcal{G}}^{k})}^{\top}\mathbf{W}^{\top}\mathbf{x}_{\mathcal{I}}^{k}+\boldsymbol{\beta}_{\mathcal{I}}^{\top}\mathbf{x}_{\mathcal{I}}^{k}+\boldsymbol{\beta}_{\mathcal{G}}^{\top}\mathbf{x}_{\mathcal{G}}^{k}+\beta_{0}\right)^{y^{k}}$
$\displaystyle\displaystyle\quad\prod_{k=1}^{N}\left[1-\sigma\left({(\mathbf{x}_{\mathcal{G}}^{k})}^{\top}\mathbf{W}^{\top}\mathbf{x}_{\mathcal{I}}^{k}+\boldsymbol{\beta}_{\mathcal{I}}^{\top}\mathbf{x}_{\mathcal{I}}^{k}+\boldsymbol{\beta}_{\mathcal{G}}^{\top}\mathbf{x}_{\mathcal{G}}^{k}+\beta_{0}\right)\right]^{1-y^{k}}$
$\displaystyle\displaystyle\quad\left(\prod_{i=1}^{|\mathcal{I}|}\prod_{\ell=1}^{L}\text{e}^{-\lambda_{W}{\theta}_{\mathcal{G}_{\ell}}\|\mathbf{W}_{i,\mathcal{G}_{\ell}}\|_{2}}\right)\times\text{e}^{-\lambda_{\mathcal{I}}\|\boldsymbol{\beta}_{\mathcal{I}}\|_{2}^{2}}\times\left(\prod_{\ell=1}^{L}\text{e}^{-\lambda_{\mathcal{G}}{\theta}_{\mathcal{G}_{\ell}}\|{\boldsymbol{\beta}}_{\mathcal{G}_{\ell}}\|_{2}}\right)$

The _maximum a posteriori_ estimation is given by:

$\displaystyle\displaystyle(\widehat{\mathbf{W}},\widehat{\boldsymbol{\beta}}_{\mathcal{I}},\widehat{\boldsymbol{\beta}}_{\mathcal{G}},\widehat{\beta}_{0})$
$\displaystyle\displaystyle\in\mathop{\text{argmax}}_{\mathbf{W},\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}},\beta_{0}}p(\mathbf{W},\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}},\beta_{0}|Y,X_{\mathcal{I}},X_{\mathcal{G}};\lambda_{W},\lambda_{\mathcal{I}},\lambda_{\mathcal{G}},\mathcal{G},\boldsymbol{\theta}_{\mathcal{G}})$
$\displaystyle\displaystyle\in\mathop{\text{argmax}}_{\mathbf{W},\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}},\beta_{0}}p(\mathbf{W},\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}},\beta_{0},Y,X_{\mathcal{I}},X_{\mathcal{G}};\lambda_{W},\lambda_{\mathcal{I}},\lambda_{\mathcal{G}},\mathcal{G},\boldsymbol{\theta}_{\mathcal{G}})$

It is equivalent to minimize the function $\displaystyle S$ defined by:

$\displaystyle\displaystyle
S(\mathbf{W},\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}},\beta_{0})$
$\displaystyle\displaystyle=-\log
p(Y,\mathbf{W},\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}},\beta_{0},X_{\mathcal{I}},X_{\mathcal{G}};\lambda_{W},\lambda_{\mathcal{I}},\lambda_{\mathcal{G}},\mathcal{G},\boldsymbol{\theta}_{\mathcal{G}})$
$\displaystyle\displaystyle=R_{N}(\mathbf{W},\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}},\beta_{0})+\Omega(\mathbf{W},\boldsymbol{\beta}_{\mathcal{I}},\boldsymbol{\beta}_{\mathcal{G}})$

## References

  * [1] Hastie, T., Tibshirani, R., Wainwright, M.: Statistical Learning with Sparsity - The Lasso and Generalizations, vol. 143. CRC Press, Boca Rato (2015)
  * [2] Ming, Y., Yi, L.: Model selection and estimation in regression with grouped variables. J. R. Statist. Soc. B, part 1 68, 49–67 (2006)
  * [3] Meier, L., van de Geer, S., Bühlmann, P.: The group lasso for logistic regression. J. R. Statist. Soc. B 70, 53–71 (2008)
  * [4] Jacob, L., Obozinski, G., Vert. J.-P.: Group lasso with overlap and graph lasso. In: Proceedings of the 26th International Conference on Machine Learning (2009)
  * [5] Beck, A., Teboulle, M.: Gradient-based algorithms with applications to signal recovery problems. In: Palomar, D.P., Eldar, Y.C. (eds.) Convex Optimization in Signal Processing and Communications, pp. 42–88. Cambribge University Press, Cambridge (2010)
  * [6] Gönen, M., Alpaydin, E.: Multiple kernel learning algorithms. J. Mach. Learn. Res. 12, 2211–2268 (2011)
  * [7] Kloft, M., Brefeld, U., Sonnenburg, S., Zien, A.: ℓp-norm multiple kernel learning. J. Mach. Learn. Res. 12, 953–997 (2011)
  * [8] Liu, J., Calhoun, V.D.: A review of multivariate analyses in imaging genetics. Front. Neuroinform. 8(29), 1–11 (2014)
  * [9] Batmanghelich, N.K., Dalca, A., Quon, G., Sabuncu, M., Golland, P.: Probabilistic modeling of imaging, genetics and diagnosis. IEEE TMI 35, 1765–1779 (2016)
  * [10] Kohannim, O., et al.: Discovery and replication of gene influences on brain structure using LASSO regression. Front. Neurosci. 6(115) (2012)
  * [11] Silver, M., et al.: Fast identification of biological pathways associated with a quantitative trait using group lasso with overlaps. Stat. Appl. Genet. Mol. Biol. 11(1), 1–40 (2012)
  * [12] Silver, M., Janousova, E., Hua, X., Thompson, P.M., Montana, G., and ADNI: Identification of gene pathways implicated in Alzheimer’s disease using longitudinal imaging phenotypes with sparse regression. NeuroImage 63(3), 1681–1694 (2012)
  * [13] Wang, H., Nie, F., Huang, H., Risacher, S.L., Saykin, A.J., Shen, L.: Identifying disease sensitive and quantitative trait-relevant biomarkers from multidimensional heterogeneous imaging genetics data via sparse multimodal multitask learning. Bioinformatics 28(12), 127–136 (2012)
  * [14] Peng, J., An, L., Zhu, X., Jin, Y., Shen, D.: Structured sparse kernel learning for imaging genetics based Alzheimer’s Disease Diagnosis. In: Ourselin, S., Joskowicz, L., Sabuncu, M.R., Unal, G., Wells, W. (eds.) MICCAI 2016. LNCS, vol. 9901, pp. 70–78. Springer, Cham (2016). doi:10.1007/978-3-319-46723-8_9
  * [15] Aiolli, F., Donini, M.: EasyMKL: a scalable multiple kernel learning algorithm. Neurocomputing 169, 215–224 (2015)
  * [16] Lorenzi, M., Gutman, B., Hibar, D., Altmann, A., Jahanshad, N., Thompson, P.M., Ourselin, S.: Partial least squares modelling for imaging-genetics in Alzheimer’s Disease: plausibility and generalization In: IEEE ISBI (2016)
  * [17] Du, L., et al.: A novel structure-aware sparse learning algorithm for brain imaging genetics. In: Golland, P., Hata, N., Barillot, C., Hornegger, J., Howe, R. (eds.) MICCAI 2014. LNCS, vol. 8675, pp. 329–336. Springer, Cham (2014). doi: 10.1007/978-3-319-10443-0_42


# Understanding Graph Convolutional Networks for Text Classification

Soyeon Caren Han,1*$\dagger$ Zihan Yuan,2* Kunze Wang,2 Siqu Long,2 Josiah
Poon1

###### Abstract

Graph Convolutional Networks (GCN) have been effective at tasks that have rich
relational structure and can preserve global structure information of a
dataset in graph embeddings. Recently, many researchers focused on examining
whether GCNs could handle different Natural Language Processing tasks,
especially text classification. While applying GCNs to text classification is
well-studied, its graph construction techniques, such as node/edge selection
and their feature representation, and the optimal GCN learning mechanism in
text classification is rather neglected. In this paper, we conduct a
comprehensive analysis of the role of node and edge embeddings in a graph and
its GCN learning techniques in text classification. Our analysis is the first
of its kind and provides useful insights into the importance of each graph
node/edge construction mechanism when applied at the GCN training/testing in
different text classification benchmarks, as well as under its semi-supervised
environment.

††* Equal contribution ††$\dagger$ Corresponding author
<EMAIL_ADDRESS>
##  1 Introduction

After the rise of deep learning, text classification models mostly applied
sequence-based learning models, CNN or RNN, which mainly captures text
features from local consecutive word sequences, but may easily ignore global
word co-occurrence in a corpus which carries non-consecutive and long-distance
semantics. Graph-based learning models are directly dealing with complex
structured data and prioritising global features exploitation. Several recent
research efforts on investigating Graph Convolutional Networks (GCN) on NLP
tasks include application to text classification (Huang et al. 2019; Yao, Mao,
and Luo 2019; Liu et al. 2020). This is largely because they can analyse rich
relational structure and preserve the global structure in graph embeddings.
The GCN-based text learning should include two main phases: 1) graph
construction from free text and 2) graph-based learning with the constructed
graph. A straightforward manner of graph construction is to represent
relationships between words/entities in the free text. Yao, Mao, and Luo
(2019) proposed a text graph-based neural network, named TextGCN, the first
corpus-level graph-based transductive text classification model. In TextGCN, a
single large textual graph is firstly constructed based on the entire corpus
with words and documents as nodes, and co-occurrence relationship between
words and documents as edges. Then, a GCN is employed to learn on the
constructed text graph. More recent studies applied extra contextual
information, such as topic model (Huang et al. 2019), syntactic and semantic
information (Liu et al. 2020), pre-trained language model (Zhang et al. 2020b)
or utilised different information propagation mechanisms (Wu et al. 2019; Zhu
and Koniusz 2021). We noticed that most studies only focus on either
hyperparameter testing or performance comparison with other state-of-the-art
text classification baselines. It is still unclear what factors in textual
graph construction or graph learning are having an impact on the GCN-based
text classification. Thus, finding the optimal textual graph construction or
learning mechanism itself, the two main phases for GCN-based text learning,
remains a black box to us. Such observations and limitations lead to several
important questions. First, the performance of GCN-based text learning methods
is highly affected by the quality of the input graph, which covers the global
structure and relations of an entire corpus or a whole dataset. Our first
question is ‘What is the best textual graph construction approach to
understand and represent the whole textual corpus?’. In a text corpus, we have
two main components, documents and words, which can be used as nodes. Then,
what feature/embedding is better to represent the node feature for the textual
graph? And what edge (relation) information should be used between nodes?
Secondly, we use GCN learning in order to capture information from the
neighbours of each word or document node. Our second question would be ‘How
much larger of a neighborhood’s information should be integrated in order to
produce the better text classification performance?’ In other words, how many
GCN layers should be stacked for the best performance on different text
classification tasks?

In this work, we focus on answering the above questions. We report the effect
of graph construction mechanisms by analysing the variants of defining main
components in a graph, including nodes and edges. Then, we present a study to
figure out the effect of GCN learning layers by integrating a variant range of
neighbourhoods’ information. We conduct our evaluation on both a full corpus
environment and a semi-supervised limited environment. The full corpus
environment is exactly the same as the original training-testing split of five
widely used benchmark corpora including 20NG, R8, R52, Ohsumed and MR, tested
by different GCN-based text classification studies (Yao, Mao, and Luo 2019;
Liu et al. 2020; Wu et al. 2019). The purpose of traditional GCN models (Kipf
and Welling 2016) is to solve semi-supervised classification tasks, so we test
on a semi-supervised environment with a very limited amount of labelled data.
For this limited setup, we use the above five widely used text classification
corpora, as well as four low resource language document classification
benchmarks (incl. Chinese, Korean, African). Note that this paper aims to
analyse the graph construction and learning mechanism of GCN-based text
classification models when there are no extra resources. This is because high-
quality resources are not always available, especially for the low resource
language or specific domain that requires expertise.

In summary, the main contributions are as follows:

  * •

We conduct a comprehensive analysis of the role of graph construction and
learning in GCN-based text classification over five widely used benchmarks
(full corpus environment) and nine benchmarks (limited training environment),
including low resource languages

  * •

We perform a comparative analysis of the accuracy performance of different
node and edge constructions when GCN is applied for text classifications

  * •

We evaluate the performance of GCN-based Text Classification with different
variants of GCN layer stacks

  * •

We make source code publicly available to encourage reproduction of the
results111https://github.com/usydnlp/TextGCN_analysis

##  2 Related Work

###  2.1 Graph Neural Networks in NLP

Graph Neural Networks (GNN) have received increasing attention in the realm of
semi-supervised learning (Kipf and Welling 2016; Li, Han, and Wu 2018).
Bastings et al. (2017) took word representations produced based on syntactic
dependency trees as graph nodes and applied them to GCN learning for machine
translation. Tu et al. (2019) proposed a Heterogeneous Document-Entity graph
and utilized a GCN to do reasoning over the constructed graph for multi-hop
reading comprehension problems. Cao et al. (2019) designed a Multi-channel
Graph Neural Network that learned the alignment-oriented knowledge graph
embeddings for entity alignment. Xu et al. (2019) extracted node features from
neighbourhoods and applied dual graph-state LSTM networks to summarize graph
local features and extracted interaction features between pairwise graphs for
entity interaction prediction. Zhang et al. (2020a) proposed automatic
sentence graph learning and incorporated it with GCN for headline generation.
Dowlagar and Mamidi (2021) combined the GCN graph modelling with multi-headed
attention for code-mixed sentiment analysis.

Some recent studies applied graph neural networks for text classification by
exploring different approaches of graph structure construction learned from
the text data. Henaff, Bruna, and LeCun (2015) and Defferrard, Bresson, and
Vandergheynst (2016) simply viewed a document as a graph node. Peng et al.
(2018) proposed a sentence-based graph in order to solve a large-scale
hierarchical text classification problem. Yao, Mao, and Luo (2019) constructed
a large textual graph with word and document nodes and edge features
represented as co-occurrence statistics, PMI/TF-IDF values. SGC (Wu et al.
2019) and S2GC (Zhu and Koniusz 2021) constructed a graph as TextGCN, but
proposed different information propagation approaches. Vashishth et al. (2019)
incorporated syntactic/semantic information for word embedding training via
GCNs. Liu et al. (2020) proposed multiple aspect graphs constructed from
external resources in terms of semantic, syntactic and sequential contextual
information, which are jointly trained. When faced with low resource text
classification problems, these approaches either do not fully explore the
latent structure within the corpus data itself as they consider only the
connections between documents, or are not applicable due to lack of external
resource. Most prior studies focused on either hyperparameter testing or
performance comparison with other state-of-the-art text classification
baselines. Distinct from these works, we examine the important factors in two
main phases of GCN-based text learning, textual graph construction and graph
learning because they have a critical impact on the GCN-based text
classification performance.

##  3 GCN-based Text Classification

We consider the task of GCN-based text classification with only single-label
classification. Figure 1 visualises the architecture of typical graph-based
text classification models based on the given corpus with no extra resources.
There are various types of GCN-based Text Classification mechanisms introduced
in the field. Two commonly used mechanisms are the corpus-level (Yao, Mao, and
Luo 2019) and document/sentence-level GCN text classification models(Huang et
al. 2019). In this work, we will focus on the former models since it captures
the global structure information of a corpus/entire dataset, whereas the
latter models consider only local-level information (from a single
sentence/document). With the former approaches, we can analyse the rich
relational structure and preserve the global structure of a graph. In this
section, we give a brief overview of GCN and TextGCN, the first corpus-level
GCN-based text classification model.

Figure 1: Text Classification with GCN. A single large graph is constructed;
words and documents appear as nodes, and co-occurrence between words and
documents. Assume that we have only three classes for the document
classification

###  3.1 Graph Convolutional Networks

GCN (Kipf and Welling 2016) is a multi-layer neural network generalized from
Convolutional Neural Networks, which directly operates on the graph-structured
data and learns representation vectors of nodes based on properties of their
neighbourhoods. Formally, a GCN graph $G$ is constructed as $G=(V,E,A)$, where
$V$ $(|V|=N)$ and $E$ represents the set of graph nodes and edges respectively
while $A\in\mathcal{R}^{N\times N}$ is the graph adjacency matrix. Based on
the constructed graph $G$, the GCN learning takes in the input matrix
$H_{0}\in\mathcal{R}^{N\times d_{0}}$ containing initial $d_{0}$-dimensional
features of the $N$ nodes in $V$ and then conducts the propagation through
layers based on the rule in equation (1), which formulates the propagation
operation from layer $l$ to the subsequent layer $(l+1)$.

$\displaystyle H^{(l+1)}=f(H^{(l)},A)=\sigma(\hat{A}H^{(l)}W^{(l)})$ (1)

Here, $\hat{A}=\tilde{D}^{-\frac{1}{2}}\tilde{A}\tilde{D}^{-\frac{1}{2}}$ is
the normalized symmetric adjacency matrix $\tilde{A}=A+I$ ($I$ refers to an
identity matrix for including self-connection of nodes); $\tilde{D}$ is the
diagonal node degree matrix, i.e. $\tilde{D}(i,i)=\sum_{j}\tilde{A}(i,j)$;
$W^{(l)}\in\mathcal{R}^{d_{l}\times d_{l+1}}$ denotes the layer-specific
trainable weight matrix for the $l$th layer ($d_{l}$/$d_{l+1}$ is the feature
dimension of layer $l$/$l+1$); $\sigma$ is a non-linear activation function
such as ReLU or softmax, which can be different for a specific layer. The main
focus of our analysis lies in exploring the role of node and edge that
constructs the graph $G$ as well as the variants of GCN learning techniques
when applying to text classification.

###  3.2 TextGCN

Inspired by a GCN, TextGCN (Yao, Mao, and Luo 2019) constructs an entire
corpus based graph, which uses all the words and documents in the corpus as
graph nodes and sets the word-word and word-document edges to preserve the
global word co-occurrence and word-document relations in the graph structure.
Then, it would be modelled by GCN learning. The edge between each word pair is
represented by the point-wise mutual information (PMI) value, normally used
for measuring semantic similarity in the Term-Sentence-Matrix. The word-
document edge is calculated based on the Term Frequency-Inverse Document
Frequency(TF-IDF) weight of the word in the document. The constructed graph is
fed into a two layer GCN as in equation (1) where the second layer node
embeddings for both word and document have the same size as the label set and
are passed into a softmax classifier for the output. The cross-entropy loss is
then calculated over all labelled documents for training and optimization.
Especially, they simply set the initial input word/document node features as
one-hot vectors.

##  4 GCN Analysis on Text Classification

###  4.1 Graph Node Construction Analysis

Following TextGCN, we set all the words and documents in the corpus as our
node set for graph $G$, i.e. the number of nodes $|V|=N=D+M$ equals to the sum
of the number of documents (corpus size $D$) and the number of unique words in
the corpus (vocab size $M$). With those word and document nodes, we explore
the role of initial node representation in GCN-based text classification
models with two commonly used input embedding types in NLP: 1) one-hot and 2)
BERT embeddings. 1) one-hot embedding is the most widely used categorical
input encoding approach in traditional NLP. 2) Bert embedding is one of the
most popular contextual word embeddings so we generate word/document node
representation. It includes each individual word embedding and a [CLS] token
for representing sentence/document-level embedding.

Figure 2: Different Variants of Edge Construction in an entire corpus-based
graph. Assume that we have four documents and three words.

$\displaystyle B_{D_{d}}=\text{BERT}(D_{d})$ (2) $\displaystyle
H_{node_{D_{d}}}=B_{D_{d}}^{[CLS]}$ (3) $\displaystyle
H_{node_{W_{m}}}=\smash{\displaystyle\min_{d\in
D_{W_{m}}}}(B_{D_{d}}^{W_{m}})$ (4)

For one-hot embedding, we simply set the feature matrix $H_{0}$ as an identity
matrix $I\in\mathcal{R}^{N\times N}$ for the one-hot vector input. For BERT
embedding, the calculation for word and document nodes are illustrated in
Equations (2)-(4). Concretely, we feed each document $D_{d}$ into the BERT
model as in equation (2), resulting in the sequence representation
$B_{D_{d}}$. For example, a document $D_{d}$ such as “John feels happy” will
result in the $B_{D_{d}}$ as “$B_{D_{d}}^{[CLS]}$ $B_{D_{d}}^{John}$
$B_{D_{d}}^{feels}$ $B_{D_{d}}^{happy}$ $B_{D_{d}}^{[SEP]}$”. We directly take
the [CLS] representation $B_{D_{d}}^{[CLS]}$ as the node embedding for
${D_{d}}$, as in equation (3). Then for a word $W_{m}$, we collect all the
documents containing this word, denoted as $D_{W_{m}}$, and apply min pooling
over all the BERT representation $B_{D_{d}}^{W_{m}}$ for this word from
documents in $D_{W_{m}}$, as is illustrated in equation (4). The essential
difference between these two types of embedding is that one-hot embedding
incorporates no external knowledge or semantic information but purely
indicates which word or document in the corpus the node is representing.
Comparatively, BERT embedding is the output representation from the BERT model
pretrained on large text corpus, which can impart the common sense semantic
information to the represented nodes and differentiate each document node
based on the document-specific word context. The detailed analysis of these
two embedding types is provided in Section 6.1.

###  4.2 Graph Edge Construction Analysis

For edge construction, we intend to fully analyse all possible co-occurring
relations between every two types of nodes, which no studies have yet
explored. We utilise document-document edges in addition to word-word and
word-document edges as can be seen in Figure 2. The construction details are
provided in equation (5). Refer to the TextGCN(Yao, Mao, and Luo 2019), we use
the PMI value between word pairs as a word-word edge feature, and utilise Term
Frequency-Inverse Document Frequency(TF-IDF) weight to weight word-document
edges.

$\displaystyle\text{PMI}(i,j)=\text{max}(\text{log}\frac{p(i,j)}{p(i)p(j)},0)$
(5) $\displaystyle p(i,j)=\frac{\text{\\#}W(i,j)}{\text{\\#}W}$ (6)
$\displaystyle p(i)=\frac{\text{\\#}W(i)}{\text{\\#}W}$ (7)

###  4.3 Graph Learning

As mentioned before, GCN only captures information about the immediate
neighbours with one layer of graph convolution. When multiple GCN hidden
layers are stacked, information about larger neighbourhoods is integrated.

We adopt this GCN propagation rule in equation (1) for modelling the
constructed graph. Specifically, we explore a various number of hidden layers
$L\in\\{1,2,3,4,5\\}$ to find the optimal range of neighbours to be
integrated. For the first $l$ ($l\in\\{1,..,L-1\\}$) layers, we apply ReLU as
activation functions as in equation (10). Here $H^{(1)}$ is the initial
feature matrix of nodes using one-hot or BERT as described in Section 4.1. We
set the last layer output for both word and document node embeddings to have
the same size as the label set and apply a softmax classifier over the output
as in equation (11). The cross-entropy loss is then calculated over all the
labelled documents as in equation (12), in which $\mathcal{Y_{D}}$ is the set
of document indices with labels available and $F$ is the output feature
dimension (equals to the number of classes). $Y$ denotes the label indicator
matrix. We provide the analysis for different numbers of hidden layers in
Section 6.3.

$\displaystyle H^{(l+1)}=\text{ReLU}(\hat{A}H^{(l)}W^{(l)})$ (8)
$\displaystyle Z=\text{softmax}(\hat{A}H^{(L)}W^{((L))})$ (9) $\displaystyle
L=-\sum_{d\in\mathcal{Y_{D}}}\sum_{f=1}^{F}Y_{df}\text{ln}Z_{df}$ (10)

##  5 Experiment setup

###  5.1 Datasets

We conduct the comprehensive analysis on both a full corpus environment and a
semi-supervised limited environment. The full environment is exactly the same
as the original training-testing split of five widely used benchmark corpora
including 20NG, R8, R52, Ohsumed and MR, followed by different GCN-based text
classification studies (Yao, Mao, and Luo 2019; Liu et al. 2020; Wu et al.
2019). The limited environment aims to cover semi-supervised text
classification with a very limited amount of labelled data. We randomly sample
1% as a training set and use the remaining 99% for testing on nine benchmarks,
including the above five benchmarks and additional four low-resource language
(incl. Chinese, Korean, African) document classification datasets222Dataset
Links: 1)http://qwone.com/˜jason/20Newsgroups/
2)3)https://www.cs.umb.edu/˜smimarog/textmining/datasets/
4)http://disi.unitn.it/moschitti/corpora.htm
5)http://www.cs.cornell.edu/people/pabo/movie-review-data/
6)https://github.com/SophonPlus/ChineseNlpCorpus/
7)https://github.com/SophonPlus/ChineseNlpCorpus/blob/master/datasets/ChnSentiCorp˙htl˙all/intro.ipynb
8)https://www.kaggle.com/captainnemo9292/korean-extremist-website-womad-hate-
speech-data 9)https://github.com/praekelt/feersum-lid-shared-task. 1)The 20NG
contains 18,846 news documents in total, which are evenly categorized into 20
classes. 2)R8 and 3)R52 (all-terms version) are subsets of the Reuters 21578
dataset. R8 has 8 topic categories with 7,674 documents, and R52 is based on
52 categories with 9,100 documents. 4)Ohsumed is collected from the MEDLINE,
which is a bibliographic database of biomedical information. Only single-label
classification task (7,400 documents) is selected. 5)MR is a binary sentiment
(positive and negative) classification dataset, which includes 10,622 short
movie review comments.

The following list shows four additional document classification benchmarks in
low-resource languages, including Chinese, Korean, and African. 6)Waimai is a
binary sentiment analysis dataset collected from a Chinese online food
ordering platform, which provides 11,987 Chinese comments about the food
delivery service. 7)ChSenti contains 7,766 Chinese documents of hotel service
comment with binary class. 8)KrHate provides 2,000 binary hate speech comments
collected from the Korean radical Anti-male online community, named Womad.
9)Xhosa is a Xhosa dataset from the NCHLT Text Corpora collected by South
African Department of Arts and Culture & Centre for Text Technology, which
contains 4,000 documents of 11 categories.

###  5.2 Implementation Details

Graph Node Setup We use one-hot and BERT embedding for the analysis. The
dimension of one-hot embedding corresponds to the number of nodes. For BERT
embeddings, “bert-base-uncased”(for English-based) and “bert-base-
multilingual-uncased”(for non-English-based) developed by Hugging Face (Wolf
et al. 2019) is used with the input dimension of 768. Graph Edge Setup In
order to construct the edge, the window size is set to 20 for PMI calculation
(word-word edges). The threshold 0.2 is applied when calculating Jaccard
similarity measure for doc-doc edges. Graph Learning Setup Each hidden layer’s
dimension is defined as 200, and the dimension of output layer is the number
of classes. The training hyperparameters include: 0.02 as the learning rate;
0.5 as the dropout rate; 0 as the $L_{2}$ loss weight; 200 as the maximum
number of epochs with early stopping of 10 epochs. Adam (Kingma and Ba 2015)
is used to train the model.

20NG R8 R52 Ohsumed MR Node feature onehot 0.8607 0.9692 0.9345 0.6824 0.7641
BERT 0.7206 0.9510 0.8440 0.4618 0.7821 Edge feature d2w only 0.8475 0.9493
0.9169 0.6667 0.7462 +w2w 0.8617 0.9693 0.9333 0.6841 0.7600 +w2w+d2d 0.8607
0.9692 0.9345 0.6824 0.7641

Table 1: Test accuracy by different node and edge construction variants on the
full environment

##  6 Discussion and Analysis

###  6.1 Effect of Node Embedding

The upper block in Table 1 shows the test accuracy by using either one-hot or
BERT embeddings as initial node features on the five benchmark datasets under
full environment. First, it can be seen that only MR achieves better result
with BERT embedding than one-hot embedding, which might be attributed to the
fact that MR as a sentiment analysis task benefits better from the general
semantics learned from a large external text. In addition, for the other four
datasets with higher accuracy via one-hot embedding, we found that datasets
with larger size of classification categories, i.e. 20NG, R52, Ohsumed, tend
to generate a bigger performance gap between the embedding types compared to
R8. We suppose that a smaller size of classification category may exert itself
better on pretrained embeddings since it does not have sufficient information
to train the global information. We also provide the comparative results under
limited environment on all the 9 datasets in Table 2 (Appendix). It is
reasonable to see an overall performance drop compared to the full setting.
For the first five benchmark datasets, MR still prefers BERT embedding while
R8 changes the preference to BERT from one-hot embedding. Both these have a
relatively small number of classification categories. This further supports
the claim that small-category datasets benefit more from BERT than large-
category counterparts. When it comes to the low resource datasets (Waimai,
ChSenti, KrHate, Xhosa), it can be seen that all of them produce better
accuracy with one-hot embedding than BERT embedding, which might be due to the
quality of pre-training on low resource language corpora. Note that we used
the edge set (d2w+w2w+d2d) for the node construction testing since it produces
the overall highest performance in both full and limited environment.

###  6.2 Effect of Edge Construction

We also analyse the usage of different edge features for both full and limited
environment in the bottom half of Table 1 and 2 respectively. More
specifically, three types of edge features are evaluated: (1) d2w only,
utilises only the word-doc edges in the constructed graph; (2) +w2w, uses both
word-doc and word-word edges; (3) +w2w+d2d, apply all the three types of edges
including doc-doc edges. Similar patterns can be found in both settings.
Overall, a d2w-only graph always results in the lowest performance, implying
insufficient global structural information conveyed by only the word-document
co-occurrence. With the w2w edge, the performance increased in all datasets
and the amount of increase mostly varies around 0.01 to 0.04 in the two
settings. In addition, d2w+w2w+d2d using all the three types of edges, further
rises the accuracy for R52 and MR under full environment and for most of the
datasets under the limited environment. This highlights the benefit of using
full co-occurring relationships (with all types of edges) in a entire graph.
The similar trend can be further observed with the persistent spatial gap
between the lowest blue line (d2w only) and the other two lines in Figure 3
(Appendix). It illustrates the corresponding accuracy for the three edge
features on the five benchmark datasets when increasing the training
proportion from extremely few labelled setting (1%) to 99%.

###  6.3 Effect of GCN Learning

The main aim of GCN learning is to capture information about immediate
neighbours with a layer of convolution. When multiple GCN layers are stacked,
information from much larger neighbourhoods is extracted and integrated. In
Figure 4 (Appendix), we conducted the text classification evaluation to find
the optimal range of neighbours’ information about each node. We stacked 1 to
5 GCN layers on different text classification in the full environment. It can
be seen that the highest performance is achieved by using 2 GCN layers for all
five datasets and the performance drops down as the layer decreases or
increases. This indicates capturing 2 levels of neighbourhood nodes is the
best and increasing the level of neighbourhoods will gradually lead to
indifferentiable node representation. 20NG and MR have a similar overall trend
and perform more consistently in the three evaluation metrics. Comparatively,
the other three datasets are observed to have much lower overall Macro F1 than
Accuracy/Weighted F1. Those trends can also be found when switching from the
full to the limited environment in Figure 5 (Appendix). Even though low
resource language datasets have extremely few labelled data, it is still 2
layers that performs the best overall, shown in Figure 5 (Appendix). When
layer number increases, the performance does not always decrease sharply as in
previous cases. Specifically, performance of the two Chinese datasets ChSenti
and Waimai goes up again at 4 layers after the decrease at 3 layers. Xhosa
only achieved a rather stable performance degradation when increasing from 2
to 5 layers. It can be seen that different languages may preserve different
patterns on the metrics with the change of GCN layers.

##  7 Conclusions

We focused on understanding the underlying factors that may influence GCN-
based text classification, and proposed to examine graph construction and
graph learning mechanisms. We systematically examined the role of node and
edge in a corpus-level textual graph, and found the optimal range of
neighbours’ information by testing the different number of GCN layer stacks.
The empirical results of experiments on various real datasets in both full
environment and limited training environment supported our analysis.

##  Appendix A Appendix

Table 2, Figure 3,4, and 5 can be found in the next page. Once the paper is
accepted, we will add those figures and tables in the main content with 6
pages (one additional page).

## References

  * Bastings et al. (2017) Bastings, J.; Titov, I.; Aziz, W.; Marcheggiani, D.; and Sima’an, K. 2017.  Graph Convolutional Encoders for Syntax-aware Neural Machine Translation.  In _Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing_ , 1957–1967.
  * Cao et al. (2019) Cao, Y.; Liu, Z.; Li, C.; Li, J.; and Chua, T.-S. 2019.  Multi-Channel Graph Neural Network for Entity Alignment.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , 1452–1461.
  * Defferrard, Bresson, and Vandergheynst (2016) Defferrard, M.; Bresson, X.; and Vandergheynst, P. 2016.  Convolutional Neural Networks on Graphs with Fast Localized Spectral Filtering.  In _NIPS_.
  * Dowlagar and Mamidi (2021) Dowlagar, S.; and Mamidi, R. 2021.  Graph convolutional networks with multi-headed attention for code-mixed sentiment analysis.  In _Proceedings of the First Workshop on Speech and Language Technologies for Dravidian Languages_ , 65–72.
  * Henaff, Bruna, and LeCun (2015) Henaff, M.; Bruna, J.; and LeCun, Y. 2015.  Deep convolutional networks on graph-structured data.  _arXiv preprint arXiv:1506.05163_.
  * Huang et al. (2019) Huang, L.; Ma, D.; Li, S.; Zhang, X.; and Houfeng, W. 2019.  Text Level Graph Neural Network for Text Classification.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , 3435–3441.
  * Kingma and Ba (2015) Kingma, D. P.; and Ba, J. 2015.  Adam: A Method for Stochastic Optimization.  _CoRR_ , abs/1412.6980.
  * Kipf and Welling (2016) Kipf, T. N.; and Welling, M. 2016.  Semi-Supervised Classification with Graph Convolutional Networks.
  * Li, Han, and Wu (2018) Li, Q.; Han, Z.; and Wu, X.-M. 2018.  Deeper insights into graph convolutional networks for semi-supervised learning.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 32.
  * Liu et al. (2020) Liu, X.; You, X.; Zhang, X.; Wu, J.; and Lv, P. 2020.  Tensor graph convolutional networks for text classification.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, 8409–8416.
  * Peng et al. (2018) Peng, H.; Li, J.; He, Y.; Liu, Y.; Bao, M.; Wang, L.; Song, Y.; and Yang, Q. 2018\.  Large-scale hierarchical text classification with recursively regularized deep graph-cnn.  In _Proceedings of the 2018 world wide web conference_ , 1063–1072.
  * Tu et al. (2019) Tu, M.; Wang, G.; Huang, J.; Tang, Y.; He, X.; and Zhou, B. 2019.  Multi-hop Reading Comprehension across Multiple Documents by Reasoning over Heterogeneous Graphs.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , 2704–2713.
  * Vashishth et al. (2019) Vashishth, S.; Bhandari, M.; Yadav, P.; Rai, P.; Bhattacharyya, C.; and Talukdar, P. 2019.  Incorporating Syntactic and Semantic Information in Word Embeddings using Graph Convolutional Networks.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , 3308–3318.
  * Wolf et al. (2019) Wolf, T.; Debut, L.; Sanh, V.; Chaumond, J.; Delangue, C.; Moi, A.; Cistac, P.; Rault, T.; Louf, R.; Funtowicz, M.; Davison, J.; Shleifer, S.; von Platen, P.; Ma, C.; Jernite, Y.; Plu, J.; Xu, C.; Scao, T. L.; Gugger, S.; Drame, M.; Lhoest, Q.; and Rush, A. M. 2019.  HuggingFace’s Transformers: State-of-the-art Natural Language Processing.  _ArXiv_ , abs/1910.03771.
  * Wu et al. (2019) Wu, F.; Souza, A.; Zhang, T.; Fifty, C.; Yu, T.; and Weinberger, K. 2019.  Simplifying graph convolutional networks.  In _International conference on machine learning_ , 6861–6871. PMLR.
  * Xu et al. (2019) Xu, N.; Wang, P.; Chen, L.; Tao, J.; and Zhao, J. 2019.  MR-GNN: Multi-Resolution and Dual Graph Neural Network for Predicting Structured Entity Interactions.  In _IJCAI_.
  * Yao, Mao, and Luo (2019) Yao, L.; Mao, C.; and Luo, Y. 2019.  Graph convolutional networks for text classification.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, 7370–7377.
  * Zhang et al. (2020a) Zhang, R.; Guo, J.; Fan, Y.; Lan, Y.; and Cheng, X. 2020a.  Structure Learning for Headline Generation.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, 9555–9562.
  * Zhang et al. (2020b) Zhang, Y.; Yu, X.; Cui, Z.; Wu, S.; Wen, Z.; and Wang, L. 2020b.  Every Document Owns Its Structure: Inductive Text Classification via Graph Neural Networks.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , 334–339.
  * Zhu and Koniusz (2021) Zhu, H.; and Koniusz, P. 2021.  Simple spectral graph convolution.  In _International Conference on Learning Representations_.

20NG R8 R52 Ohsumed MR Waimai ChSenti KrHate Xhosa Node feature onehot 0.6937
0.8973 0.7990 0.4092 0.6178 0.8301 0.7548 0.9048 0.9952 BERT 0.6585 0.9147
0.7962 0.3893 0.7255 0.8059 0.6970 0.8256 0.7878 Edge feature d2w only 0.6239
0.8698 0.7776 0.3906 0.5995 0.8149 0.5297 0.7445 0.9874 +w2w 0.6680 0.8949
0.7978 0.4099 0.6158 0.8283 0.7410 0.7609 0.9953 +w2w+d2d 0.6937 0.8973 0.7990
0.4092 0.6178 0.8301 0.7548 0.9048 0.9952

Table 2: Test accuracy by different node and edge construction variants on the
limited environment Figure 3: Test accuracy by varying training data
proportions (from 1% to 99%) Figure 4: Test performance by varying GCN hidden
layer stacks on the full environment Figure 5: Test performance by varying GCN
hidden layer stacks on the limited environment