International Symposium on New Developments of Theories andMethodologies for Large Complex Data
November 5-6, 2021
Venue:
Conference Room 102, Tsukuba International Congress Center
2-20-3 Takezono, Tsukuba, Ibaraki 305-0032, Japan (Hybrid Symposium with Zoom)
Organizers:
Makoto Aoshima (University of Tsukuba)
Mika Sato-Ilic(University of Tsukuba)

Kazuyoshi Yata

(University of Tsukuba)

Aki Ishii(Tokyo University of Science)

Yugo Nakayama(Kyoto University)

Supported by
Grant-in-Aid for Scientific Research (A) 20H00576 (Project Period: 2020-2024)“Innovative Developments of Theories and Methodologies for Large Complex Data”(Principal Investigator: Makoto Aoshima)
Grant-in-Aid for Challenging Research (Exploratory) 19K22837 (Project Period: 2019-2021)“Tackling individualized modeling with ultra-high dimensional data”(Principal Investigator: Makoto Aoshima)

Program (UTC+9)
November 5 (Friday)
14:00∼14:10

Opening

14:10∼15:00 Koji Tsukuda∗,a and Shun Matsuurabab(Faculty of Mathematics, Kyushu University)(Faculty of Science and Technology, Keio University)

On high-dimensional testing for common principal components model
15:00∼15:50 Takahiro Nishiyama∗,a , Masashi Hyodob and Tatjana Pavlenkoc(Zoom)a(Department of Business Administration, Senshu University)b(Department of Economics, Kanagawa University)c(Department of Mathematics, KTH Royal Institute of Technology)

A two sample Behrens-Fisher problem for factor models in high dimensions

15:50∼16:40 Shao-Hsuan Wang
(Zoom)(Graduate Institute of Statistics, National Central University)

On the asymptotic result of Kronecker envelope principal component analysis
in high dimension low sample size data

16:45∼19:00

Special Invited and Keynote Sessions 1(Zoom)

November 6 (Saturday)
10:00∼10:50 Yugo Nakayama∗,a , Kazuyoshi Yatab and Makoto Aoshimabab(Graduate School of Informatics, Kyoto University)(Institute of Mathematics, University of Tsukuba)

Asymptotic properties of high-dimensional kernel PCA and its applications
10:50∼11:40 Thong Pham∗,a , Paul Sheridanb and Hidetoshi Shimodairac(Zoom)a(RIKEN AIP)b

(Tupac Bio)c

(Graduate School of Informatics, Kyoto University, and RIKEN AIP)

Estimating preferential attachment in growing networks
11:40∼13:00

Lunch

13:00∼13:50 Tsubasa Itoa and Shonosuke Sugasawa∗,b(Zoom)

a(M&D Data Science Center, Tokyo Medical and Dental University)b(Center for Spatial Information Science, The University of Tokyo)

Grouped generalized estimating equations for heterogeneous longitudinal data
13:50∼14:40 Shota Katayama(Zoom)(Faculty of Economics, Keio University)

Hypothesis testing on high dimensional parameter under confounding
14:40∼15:30 Tsutomu T. Takeuchi(Division of Particle and Astrophysical Science, Nagoya University)

New quantification of galaxy evolution via manifold learning

15:35∼17:50

Special Invited and Keynote Sessions 2(Zoom)
17:50∼18:00

Closing(∗ Speaker)

Special Invited and Keynote Sessions 1
November 5 (Friday)
Special Invited Talk
16:45∼17:45

Double data piling and negatively ridged classifiers in high dimensions

Speaker: Sungkyu Jung
(Department of Statistics, Seoul National University)
Chair: Kazuyoshi Yata (Institute of Mathematics, University of Tsukuba)

Keynote Talk
17:50∼19:00

Using the classical Growth Curve model in non-standard situations

Speaker: Dietrich von Rosen(Department of Energy and Technology, Swedish University of Agricultural Sciences)
Discussion Leader: Shinpei Imori(Graduate School of Advanced Science and Engineering, Hiroshima University)

Special Invited and Keynote Sessions 2
November 6 (Saturday)
Special Invited Talk
15:35∼16:35

Normal-reference tests for high-dimensional hypothesis testing

Speaker: Jin-Ting Zhang(Department of Statistics and Data Science, National University of Singapore)
Chair: Shota Katayama (Faculty of Economics, Keio University)

Keynote Talk
16:40∼17:50

Generalized information criterion for high-dimensional PCA rank selection

Speaker: Su-Yun Huang(Institute of Statistical Science, Academia Sinica)
Discussion Leader: Yuan-Tsung Chang(Department of Social Information, Mejiro University)

On high-dimensional testing for common principal components model *Koji Tsukuda andShun Matsuura

In this presentation, a high-dimensional statistical test procedure for the common principal components (CPC) hypothesis on two population covariance matrices was discussed. The CPC hypothesis meansthat the population covariance matrices Σx and Σy can be simultaneously diagonalized; that is to say,there exists an orthogonal matrix U such that both U ⊤ Σx U and U ⊤ Σy U are diagonal. Testing the CPC
hypothesis was ﬁrstly considered by Flury [2], and has been studied in several works; see, for example,
Boente, Pires, and Rodrigues [1] and Hallin, Paindaveine and Verdebout [3, 4]. Unlike these studies, weconsidered the case where sample sizes are smaller than the number of observed variables.
Our approach is based on the asymptotic normality ofM=1√
p2 n a n b n c n dtr(Ta (na )Tb (nb )Tc (nc )Td (nd )),that is the trace of the product of four p × p independent Wishart matrices Ta (na ), Tb (nb ), Tc (nc ), Td (nd )with respective degrees-of-freedom na , nb , nc , nd and respective scale matrices Σa , Σb , Σc , Σd , under ahigh-dimensional asymptotic regimena , nb , nc , nd ≍ pδ , 0 < δ < 1.
The test statistic is shown to be asymptotically distributed as the standard normal distribution under
the null hypothesis. It is also shown that the asymptotic power of the test goes to 1 under the alternative
hypothesis. The performances of the test for ﬁnite samples are numerically examined.
This presentation was based on the article [5].

*This work is supported by JSPS KAKENHI (grant number: 18K13454,20K11713,21K13836)
Faculty of Mathematics, Kyushu University
Faculty of Science and Technology, Keio University

1References[1] Boente, G., Pires, A.M., and Rodrigues, I.M. (2009). Robust tests for the common principal components model. J. Statist. Plann. Inference 139, no.4, 1332–1347.
[2] Flury, B.N. (1984). Common principal components in k groups. J. Amer. Statist. Assoc. 79, no.388,
892–898.
[3] Hallin, M., Paindaveine, D., Verdebout, T. (2010). Testing for common principal components underheterokurticity. J. Nonparametr. Stat. 22, no.7, 879–895.
[4] Hallin, M., Paindaveine, D., Verdebout, T. (2013). Optimal rank-based tests for common principal
components. Bernoulli 19, no.5B, 2524–2556.
[5] Tsukuda, K., Matsuura, S. (2021). Limit theorem associated with Wishart matrices with application
to hypothesis testing for common principal components. J. Multivariate Anal. 186, 104822.

2

A two sample Behrens-Fisher problem for factor models in highdimensions
Takahiro Nishiyamaa , Masashi Hyodob and Tatjana PavlenkocaDepartment of Business Administration, Senshu Universityb
Faculty of Economics, Kanagawa Universityc
Department of Mathematics, KTH Royal Institute of Technology
Let xgi “ pxgi1 , . . . , xgip qJ „ Fg be iid p-dimensional random vectors collected fromthe ith subject in the gth population, where Fg denotes the distribution function for gth
population, i P t1, . . . , ng u, g P t1, 2u. A factor model assumes that for each g P t1, 2u, theobservable vector xgi is decomposable into a latent factor and an idiosyncratic component
as follows:
xgi “ µg ` Bg zgi ` Ψ1{2g ϵgi ,

(1)

where µg P Rp is a deterministic intercept vector, zgi “ pzgi1 , . . . , zgidg qJ is the dg dimensional latent factor vector, and ϵgi “ pϵgi1 , . . . , ϵgip qJ is the p-dimensional error
vector which is uncorrelated with the latent factor. In what follows, we assume thatdg P N is a fixed number. Further, Bg “ pbg1 , . . . , bgp qJ denotes the loading matrix
where for each j P t1, . . . , pu, bgj “ pbgj1 , . . . , bgjdg qJ P Rdg is a non-random vector, andΨg “ diagpψg1 , . . . , ψgp q is the non-random p ˆ p diagonal matrix whose elements areψg1 ą 0, . . . , ψgp ą 0. For the latent vector zgi and error vector ϵgi , we further assume
that zgiℓ are iid with Epzgiℓ q “ 0, Epz2giℓ q “ 1 and Epz4giℓ q “ κzg ă 8, and ϵgij are iid with
Epϵgij q “ 0, Epϵ2gij q “ 1 and Epϵ4gij q “ κϵg ă 8 for g P t1, 2u, i P t1, . . . , ng u, j P t1, . . . , pu
and ℓ P t1, . . . , dg u. Structural assumptions of the model (1) imply that
Epxgi q “ µg , covpxgi q “ Bg BJ
g ` Ψg :“ Σg ,pˆpwhere Σg P Rpˆp
ą0 and Rą0 denotes the space of real, symmetric, positive definite, p ˆ pmatrices.
By using the data generated by (1), we designed a high-dimensinal test procedure fortesting

H : µ1 “ µ2 , A : µ1 ‰ µ2 .

(2)The proposed test will be valid for high-dimensional, non-normal, unbalanced data undertwo-sample Behrens-Fisher setting for the low-dimensional factor model.
Recently, a number of useful two-sample tests have been proposed for high-dimensionalsettings. Nevertheless, their applicability relies on a critical assumption, that is, requiring
elements of xgi to be weakly dependent. Being expressed in terms of the latent factor
model (1) and assuming the common covariance matrix Σ, the weak dependence is equivalent to an assumption stated as trpΣ4 q{tr2 pΣ2 q Ñ 0, when p Ñ 8; see e.g. assumption

(3.6) in Chen and Qin (2010). This assumption is crucial for establishing the asymptotic
normality of the test statistic proposed by Chen and Qin (2010); see Theorem 1 of theirpaper. However, besides of being diﬃcult to verify in practice, the regularity assumptions
imposed on the covariance structure in e.g. Bai and Saranadasa (1996) and Chen andQin (2010) can be easily violated in covariance models where the eigenvalues of Σg are
dominated by few top ones. This is precisely the type of covariance structure underlying
xgi : under the low-dimensional latent factor model (1), the first dg eigenvalues of Σg are
considerably larger than than the rest. So, Ma et al (2015) proposed the mean diﬀerence
test for the model (1) under homoscedastisity. However, an assumption of common covariance matrix is a very strong assumption which is hard to be practically verified in p " ngsettings, this in turn limits the applicability of the procedure developed by these authors.
We do not assume that Σ1 “ Σ2 ; out test statistics, along with their limit properties arestudied under heteroscedasticity, i.e. solves a general, two-sample Behrens-Fisher problem
for the latent factor model (1).
For this problem, we defined the data-driven test statistic as
"
*
n
1 {
1 {
2
TF A “}x1 ´ x2 } ´ trpΨ1 q ´ trpΨ2 q ,
p
n1
n2
řdpg
řng
{
xgi and trpΨ
where, for g P t1, 2u, xg “ p1{ng q i“1g q “ trpSg q ´ℓ“1 λℓ pSg q. Here,řng
λℓ pSg q is the ℓth largest eigenvalue of matrix Sg “ t1{png ´ 1qu i“1 pxgi ´ xg qpxgi ´ xg qJand dpg is a consistent estimator of dg based on the ER method proposed by Ahn andHorenstein (2013). Besides, we derived the limiting null distribution of TF A under someassumptions and constructed test procedure for testing p2q. Also, we compared, throughsimulations, the performance of the proposed test and existing procedures suitable for a
two-sample, Behrens-Fisher problem in high-dimensional data in terms of size control andpower.
References[1] Ahn, S. C., Horenstein, A R., 2013. Eigenvalue ratio test for the number of factors. Econometrica, 81, 1203–1227.
[2] Bai, Z.D., Saranadasa, H., 1996. Eﬀect of high dimension: by an example of a two sample
problem. Statist. Sinica, 6, 311–329.
[3] Chen, S.X., Qin, Y.L., 2010. A two-sample test for high dimensional data with applications
to gene-set testing. Ann. Statist., 38, 808–835.
[4] Ma, Y., Lan, W., Wang, H., 2015. A high dimensional two-sample test under a low dimensional factor structure. J. Multivariate Anal., 140, 162–170.

On the asymptotic result of Kronecker envelope
principal component analysis in high dimension low
sample size data
Shao-Hsuan Wang
(Graduate Institute of Statistics, National Central University)

Abstract
In the analysis of high dimensional data, the Kronecker envelope principal component
analysis (KEPCA) serves as an eﬃcient alternative to the classical PCA at pre-processing
steps. We derive the consistency and the asymptotic normality of Kronecker envelope
PCA in High dimension low sample size (HDLSS) setting and compare it with classical
PCA. Face database example is used to demonstrate our method.

1

Double data piling and negatively ridged classifiersin high dimensions
Sungkyu Jung
(Department of Statistics, Seoul National University)

Abstract
Data piling refers to the phenomenon that training data vectors from each class project
to a single point for classification. While this interesting phenomenon has been a key tounderstanding many distinctive properties of high-dimensional discrimination, the theoretical underpinning of data piling is far from properly established. In this work, highdimensional asymptotics of data piling is investigated under a spiked covariance model,
which reveals its close connection to the well-known ridged linear classifier. In particular,
by projecting the ridge discriminant vector onto the subspace spanned by the leading
principal component directions and the maximal data piling vector, we show that a negatively ridged discriminant vector can asymptotically achieve data piling of independent
test data, essentially yielding a perfect classification. The second data piling direction is
obtained purely from training data and shown to have a maximal property. Furthermore,
asymptotic perfect classification occurs only along the second data piling direction. This
interesting phenomenon is shown to also occur in multi-category classification problems,
in which the second data piling subspaces are estimated by negatively ridged discriminant
subspaces. We demonstrate that negative ridge parameters can be optimal in classificationof well-known image and microarray datasets.

1

Using the classical Growth Curve model in
non-standard situations
Dietrich von Rosen(Department of Energy and Technology, Swedish University of Agricultural Sciences)

Abstract
More than ﬁfty years ago the Growth Curve model (GMANOVA) was introduced. It
is a multivariate bilinear regression model which under a normality assumption belongs
to the curved exponential family. The model is of great interest since explicit maximum
likelihood estimators exist. In the presentation we will focus on estimation based on thelikelihood. Four diﬀerent types of models will be discussed via the Growth Curve modeland all types are connected to high-dimensional problems.
The ﬁrst case is a discussion of the Growth Curve model when the number of independent observations are less than the number of repeated measurements. It turns outthat the Moore-Penrose generalized inverse of a singular Wishart matrix can be used.
Unfortunately, moments for the Moore-Penrose inverse are not available if the dispersion
matrix is unstructured and we discuss what to do in this situation.
The second case which will be considered is about estimation in Partial Least Squares(PLS). PLS is often connected with diﬀerent algorithms. In this talk, when estimating
parameters, we will set up a model which gives the same result as algorithms. In PLS a
Krylov space plays a central role.
The third case which will be discussed is when in the Growth Curve model there arerank restrictions on the mean parameters of the model. When there are a huge numberof background variables this type of model can be useful. If one can imagine that there
are a few number of latent processes which govern the background variables the rank
restrictions make sense.
In the fourth case we extend the third one by also allowing for rank restrictions on
the dispersion matrix. Rank restrictions on an unknown dispersion matrix has not beenstudied much and we will show how the mathematics for ﬁnding estimators can be carriedout. However, there are still many open problems connected to the interpretability of theresults and connected to the decision of the rank on both the mean parameters and thedispersion of the model.

1Asymptotic properties of high-dimensional
kernel PCA and its applications
Yugo Nakayamaa , Kazuyoshi Yatab and Makoto Aoshimaba

Graduate School of Informatics, Kyoto University
Institute of Mathematics, University of Tsukubab1

Introduction

In this talk, we considered the kernel principal component analysis (KPCA) for high-dimension,
low-sample-size (HDLSS) data. We proposed clustering by using asymptotic properties of the
KPCA and applied it to outlier detection.

2

Asymptotic properties of high-dimensional kernel PCA

Suppose there are independent and d-variate populations, Πi , i = 1, 2, having an unknown
mean vector µi and unknown covariance matrix Σi for each i. Suppose we have a d × n data
matrix X = (x1 , ..., xn ), where xj s are independently taken from Π1 or Π2 . Let{}
ni = # j|xj ∈ Πi for j = 1, ..., n ,where #A denotes the number of elements in a set A. Note that n = n1 + n2 . We assume
that n and ni s are independent of d, and ni ≥ 1 for i = 1, 2. For the sake of simplicity, we
assume that tr(Σ1 ) ≤ tr(Σ2 ) andxj ∈ Π1 , j = 1, ..., n1 , xj ∈ Π2 , j = n1 + 1, ..., n.

(1)

In this section, we consider asymptotic properties of the KPCA with the linear kernelk(xj , xj ′ ) = x⊤j xj ′ in the HDLSS context that d → ∞ while n is fixed. Let K be an n × ngram matrix with the (j, j ′ ) element k(xj , xj ′ ). Let P n = I n − n−1 1n 1⊤
n , where I n denotes
the n-square identity matrix and 1n = (1, ..., 1)⊤ . We define the (centroid) gram matrix by
K 0 = P n KP n .
Note that rank(K 0 ) ≤ n − 1. Let λ̂1 ≥ · · · ≥ λ̂n−1 be the eigenvalues of K 0 . Then, we define
the eigen-decomposition of K 0 by
K0 =

n−1
∑
i=1

1

λ̂i ûi û⊤
i ,where ûi = (ûi1 , ..., ûin )⊤ denotes a unit eigenvector corresponding to the λ̂i . The ith (nor√malized) PC score of xj is given by sij = nûij . Since the sign of an eigenvector is arbitrary,⊤
we assume that (1⊤n1 , −1n2 )û1 ≥ 0 without loss of generality.
We assume the following condition:(A-i) Var(∥x − µi ∥2 |x ∈ Πi ) = O{tr(Σ2i )} as d → ∞ for i = 1, 2.
Note that E(∥x − µi ∥2 |x ∈ Πi ) = tr(Σi ) for i = 1, 2. If Πi s are Gaussian, it holds that
Var(∥x − µi ∥2 |x ∈ Πi ) = 2tr(Σ2i ) for i = 1, 2, so that (A-i) naturally holds. Let ∆µ =∥µ1 − µ2 ∥2 and ∆Σ = |tr(Σ1 ) − tr(Σ2 )|, where ∥ · ∥ denotes the Euclidean norm. Here, we
assume the following conditions:(A-ii) tr(Σ2i )/∆2µ = o(1) as d → ∞ for i = 1, 2.
(A-iii) lim supd→∞∆Σ< 1 when n2 ≥ 2.
n1 ∆µNakayama et al. (2021) gave the following result.
Theorem 2.1. Assume (A-i) to (A-iii). Then, it holds that as d → ∞{√n /n + oP (1),j = 1, ..., n1 ,
√2 1
s1j =
− n1 /n2 + oP (1), j = n1 + 1, ..., n.

(2)

Thus from Theorem 2.1, one can classify xj s into two groups by the sign of the first PC
scores.
Remark 1. If ∆µ /∆Σ is small, we do not recommend to use the linear kernel. In such case,
you can use the Gaussian kernel k(xj , xj ′ ) = exp(−∥xj − xj ′ ∥2 /γ), where γ > 0. Note that
the Gaussian kernel can draw information about heteroscedasticity via the diﬀerence of Σi s.
We considered applying the proposed clustering method to outlier detection. We set
n1 = 1 and x1 is regarded as an outlier. From (2), under some regularity conditions, it holds
as d → ∞{√n − 1 + oP (1),j = 1,
s1j =√−1/ n − 1 + oP (1), j = 2, ..., n.
Thus by using the first PC scores, one can detect the outlier.

References[1] Y. Nakayama, K. Yata, M. Aoshima, Clustering by principal component analysis with
Gaussian kernel in high-dimension, low-sample-size settings, J. Multivariate Anal. 185(2021) in print.

2

Estimating preferential attachment in growing networks
Thong Pham∗1 , Paul Sheridan2 , and Hidetoshi Shimodaira3,1
1

RIKEN Center for AIP
2
Tupac Bio, Inc.
3
Kyoto University

We consider the problem of estimating the preferential attachment (PA) function Ak fromempirical data that does not contain any information about the growth process of the network.
Consider a growing network that grows from time-step t = 1 to T and denote its snapshot at timestep t by Gt . Traditionally, the estimation of Ak is often considered when the growth process of the
network can be observed at at least two time-steps [6, 5, 7]. However, what if we cannot observeanything about the growth process and have to content ourselves with only the one snapshot
GT ? When it comes to estimating PA in general growing networks without time-resolved data,no satisfactory methods exist. All existing methods assume either unrealistic network types or
unnecessarily restrictive functional forms for Ak [1, 4, 2, 3].
We propose a method called PAFit-oneshot to nonparamterically estimate the PA function of a
growing network from its final snapshot, GT , alone [9]. Our method does not assume any functional
form for the PA function, and can be applied to any real-world network snapshot.
Our method can be summarized as follows. From Theorem 1 in [9], one gets the basic equation:
P
E j>k nj (T )
,
(1)
Ak ≈
Enk (T )P
where nk (t)Pis the number of degree k nodes at time-step t. If one estimates E j>k nj (T ) andEnk (T ) by j>k nj (T ) and nk (T ), one arrives at the baseline estimator:
PÂbaseline=kj>k nj (T )nk (T ),

(2)

which is the estimator of Gao et al. [4]. However, there is a severe underestimation of Ak in theregion of large k by this method (see Fig. 1 of [9]). We call this underestimation the waterfall
artefact.
The root of the waterfall artefact is that nk (T ) is a poor estimator of Enk (T ) when k is large.
Given GT , we only estimate the Ak values for the degrees k observed in GT , which means nk (T )is positive a priori. Therefore, nk (T ) is actually an estimator for the conditional expectation
E [nk (T ) | nk (T ) > 0], which is equal to Enk (T )/P(nk (T ) > 0). Surprisingly, the presence of thisbias, let alone a proposed correction for it, has never been discussed in the literature.
∗ Corresponding author. Email: thong.pham@riken.jp

1

To remove this bias, we need to estimate pk = P(nk (T ) > 0), the probability that k exists
in GT . This resembles the problem of selective inference. Selective inference typically considers
model selection in a regression setting and adjusts the bias of regression coefficients for the selected
predictors; one must correct for the effect of choosing the predictor. In our problem, model selection
is equivalent to choosing over which values of k to use to estimate Ak , which is where k exists in
GT . Correcting for pk leads us to the following equation:
Pj>k nj (T )
.
(3)
Ak ≈
nk (T )pk
Starting from initial rough estimations of pk and Ak , our method iteratively improves theseestimations using Monte Carlo simulations (see Fig. 6 of [9]). The proposed method is implemented
in the R package PAFit [8].

References[1] Ivona Bezáková, Adam Kalai, and Rahul Santhanam. Graph model selection using maximum
likelihood. In Proceedings of the 23rd International Conference on Machine Learning, ICML
’06, pages 105–112, New York, NY, USA, 2006. ACM.
[2] George T. Cantwell, Guillaume St-Onge, and Jean-Gabriel Young. Inference, model selection,and the combinatorics of growing trees. Phys. Rev. Lett., 126:038301, Jan 2021.
[3] Harry Crane and Min Xu. Inference on the history of a randomly growing tree. Journal of theRoyal Statistical Society: Series B (Statistical Methodology), 83(4):639–668, 2021.
[4] Fengnan Gao, Aad van der Vaart, Rui Castro, and Remco van der Hofstad. Consistent estimation in general sublinear preferential attachment trees. Electron. J. Statist., 11(2):3979–3999,
2017.
[5] H Jeong, Z Néda, and AL Barabási. Measuring preferential attachment in evolving networks.
Europhysics Letters, 61(61):567–572, 2003.
[6] MEJ Newman. Clustering and preferential attachment in growing networks. Physical Review
E, 64(2):025102, 2001.
[7] Thong Pham, Paul Sheridan, and Hidetoshi Shimodaira. PAFit: A statistical method formeasuring preferential attachment in temporal complex networks. PLOS ONE, (9):e0137796, 9
2015.
[8] Thong Pham, Paul Sheridan, and Hidetoshi Shimodaira. PAFit: An R package for the nonparametric estimation of preferential attachment and node fitness in temporal complex networks.
Journal of Statistical Software, Articles, 92(3):1–30, 2020.
[9] Thong Pham, Paul Sheridan, and Hidetoshi Shimodaira. Non-parametric estimation of thepreferential attachment function from one network snapshot. Journal of Complex Networks,
9(5), 09 2021. cnab024.

2

Grouped generalized estimating equations forheterogeneous longitudinal data
Tsubasa Ito1 and Shonosuke Sugasawa2
1 M&D Data Science Center, Tokyo Medical and Dental University
2 Center for Spatial Information Science, The University of Tokyo

Longitudinal data in which responses (repeated measurements) within the same
subject are correlated is appeared in many scientiﬁc applications such as biomedical
statistics and social science. For analyzing longitudinal data, it is typically diﬃcult
to correctly specify the underlying correlation structures among response variables
within the same subject, and one of the standard approaches is the generalized estimating equations (GEE) developed by Liang and Zeger (1986), which uses “working”
correlation structures speciﬁed by users. The advantage of the GEE approach is thatthe estimator is still consistent even when the working correlation is misspeciﬁed.
However, the existing GEE methods assumes homogeneous regression coeﬃcientsthat are common to all the subject, which could be restrictive in practical applications since there might be potential heterogeneity among subjects or clusters as
conﬁrmed in several applications.
In this work, we extend the standard GEE analysis to take account of potential
heterogeneity in longitudinal data. Speciﬁcally, we develop grouped GEE analysis
by adopting the grouping approach that is widely adopted in literatures for panel
data analysis. We assume that subjects in longitudinal data can be classiﬁed to aﬁnite number of groups, and subjects within the same group share the same regression coeﬃcients, that is, the regression coeﬃcients are homogeneous over subjects inthe same groups. Since the grouping assignment of subjects are unknown, we treat
it as unknown parameters and estimate them and group-wise regression coeﬃcientsimultaneously. Given the grouping parameters, the standard GEE can be performed
to obtain group-wise estimators of regression coeﬃcients. On the other hand, given
the group-wise regression coeﬃcients, we consider estimating the grouping parameters using a kind of Mahalanobis distance between response variables and predictors
with taking account of potential correlations via working correlation matrix. In otherwords, we employ the working correlation not only in performing GEE analysis ineach group but also estimating the grouping assignment. We will show that thegrouped GEE method can be easily carried out by a simple iterative algorithm similar to k-means algorithm that combines the existing algorithm for the standard GEE
and simple optimization steps for grouping assignment. Moreover, we adopt the crossvalidation to carry out data-dependent selection of the number of groups.
We derive the statistical properties of the grouped GEE estimator in an asymptotic framework where both n (the number of subjects) and T (the number of repeatedmeasurements) tend to inﬁnity, but we here allow T to grow considerably slower thann, namely, n/T ν → 0 for some large ν. Hence, our method can be applicable when T
is much smaller n as observed in many applications using longitudinal data. As theoretical diﬃculties of the grouped estimation in longitudinal data analysis, the truecorrelations within the same subject can be considerably high, so the existing theoretical argument assuming negligibly small correlations imposed typically by mixing
conditions for the underlying true correlations is no more applicable. To overcome
the limitation of the existing theoretical argument, we consider grouping assignment
using a kind of Mahalanobis distance with working correlation, and we will showthat such grouping strategy leads to consistent estimation of the grouping parameters as long as the working correlation is relatively close to the true one. Therefore,
even when the underlying correlations within the same subject is not weak, we cansuccessfully estimate the grouping parameters using a reasonable working correlation
matrix. Then, we will also establish consistency and asymptotic normality of thegrouped GEE estimator of the regression coeﬃcients, and also provide a consistentestimator of asymptotic variances.

2

Hypothesis testing on high dimensional parameter underconfounding
Shota Katayama
Faculty of Economics, Keio University, Japan

1

Introduction

In high dimensional data analysis, especially in diﬀerential gene expression analysis, detecting
diﬀerence of a huge number of features between two groups is essential problem. Eachindividual is assigned to one of two groups according to diﬀerence in conditions, e.g., whether
they have been treated or not. In situations where one can assume a random assignment, thetwo groups are diﬀerent due to the condition alone, and thus a pure comparison is possible.
As in Potter (2003) and Heller et al. (2009), however, the random assignment is not feasible
in gene data analysis where the condition may be cancer or non-cancer. In this case, thediﬀerence between the two groups may stem from variables other than the condition ofinterest, and hence such confounding variables need to be adjusted appropriately. Moreover,gene data are not necessary continuous values, but can be binary or count.
This talk provides an uniﬁed inference on high dimensional parameter for the comparisonof two groups. Suppose we have p dimensional response vector Y = (Y1 , . . . , Yp )T , q dimensional confounding vector X = (X1 , . . . , Xq )T with X1 = 1 for intercept, and an indicatorvariable D ∈ {0, 1} introducing the diﬀerence in conditions. Assume that()
E(Yj |D, X) = h Dτj∗ + X T θj∗ ,j = 1, 2, . . . , p,P(D = 1 |X) = 1/(1 + exp(−X T γ ∗ )),

(1)
(2)

where τj∗ ∈ R, θj∗ ∈ Rq , γ ∗ ∈ Rq and h(·) is a known diﬀerentiable inverse link function. The
ﬁrst equation (1) models the dependence between Yj and (D, X) under a general framework
to deal with various type of responses. For instance, h(·) would be identity function when
Yj is continuous and would be exp(·) when Yj is count. The second (2) assumes the logistic
model for D. Finally, suppose we have n i.i.d. observations {(Yi , Di , Xi )}ni=1 following (1)and (2) and consider the high dimensional setting where p is greater than n, but q is smaller.

In (1) the parameter τj∗ represents the eﬀect of D on Yj , and thus the high dimensionalparameter vector τ ∗ = (τ1∗ , . . . , τp∗ )T is of interest. Under the strong ignorability, τ ∗ can also
be regarded as a causal parameter. This talk provides an inference on τ ∗ under the highdimensional asymptotic framework (n, p) → ∞.

2

Proposed methodology

Applying the maximum likelihood estimation to each Yj as the generalized linear regression,
we can obtain an estimator of τj∗ . In the case where X is confounded, however, the MLE ofτj∗ is not eﬃcient since it ignores the model (2). Chernozhukov et al. (2018) and Vansteelandt
and Dukes (2020) provide an eﬃcient score function of τj∗ given by
Si (τj ; θj∗ , γ ∗ ) =Yij − h(Di τj + XiT θj∗ )(Di − πi (γ ∗ )),h′ (Di τj + XiT θj∗ )(3)

where πi (γ ∗ ) = 1/(1 + exp(−XiT γ ∗ )). The nuisance parameters (θj∗ , γ ∗ ) can be estimatedby MLE or other methods and plugged into the score function above, but we need the crossﬁtting approach (Chernozhukov et al. (2018)) to control the asymptotic behavior of (3) withestimated nuisance parameters.
Based on the estimate τ̂j and its variance estimate σ̂j2 by the cross-ﬁtting, this talk gives
a maximum type test statistics for the global hypothesis H0 : τ ∗ = 0. A multiple testingprocedure with controlling the false discovery rate is also provided.

References
Chernozhukov, V., Chetverikov, D., Demirer, M., Duﬂo, E., Hansen, C., Newey, W. andRobins, J. (2018). Double/debiased machine learning for treatment and structural parameters, The Econometrics Journal, 21(1), C1–C68.
Heller, R., Manduchi, E. and Small, D. S. (2009). Matching methods for observational microarray studies. Bioinformatics, 25(7), 904–909.
Potter, J. (2003). Epidemiology, cancer genetics and microarrays: making correct inferences,using appropriate designs. Trends in Genetics, 19, 690–695.
Vansteelandt, S. and Dukes, O. (2020). Assumption-lean inference for generalised linear
model parameters. http://arxiv.org/abs/2006.08402

New Quantification of Galaxy Evolution via Manifold Learning
Tsutomu T. TAKEUCHI1,2, Suchetha COORAY1,†, Kai T. KONO1
1. Division of Particle and Astrophysical Science, Nagoya University, Nagoya 464-8602, Japan
2. The Research Center for Statistical Machine Learning, the Institute of Statistical Mathematics, Tachikawa, Tokyo 190-8562, Japan† JSPS Research Fellow (DC1)

1. Galaxy Formation and Evolution
Matter in the early Universe was almost uniform, and a slightly dense region grew by gravity,
finally into a galaxy. It was attempted to develop a theory to deal with the star formation andassociated history of heavy element synthesis, under an assumption that a galaxy has formed from asingle, huge gas cloud. While the research in this direction was once completed in the first half of
1980s, this was not the end of the studies of galaxy evolution. Cosmological research that has
progressed in parallel has revealed that galaxies merge and grow. This indicates that the galaxy
evolution is a very complicated process that strongly depends on the density of the surrounding
galaxies and the gas density. In order to formulate the galaxy evolution, it is necessary to determine
such a huge system of equations. Though astrophysicists have constructed the governing equations
from the physical laws from the first principle before, such a method is not realistic anymore whenthe quantity space exceeds 10 dimensions. Galaxy surveys as of the 21st century provides hundredsof physical quantities for hundreds of millions of galaxies, typical big data in both quality and quantity
indeed. The feature space of the galaxy to be analyzed exceeds 100 dimensions. Therefore, thecharacterization of the galaxy evolution is no longer possible by the traditional method relying on
physical intuition.

2. Galaxy Manifold
2.1 Rise, fall, and revival of the galaxy manifold
From 1970s to the mid-1980s, classical multivariate analysis methods such as the principal
component analysis (PCA) were used to combine physical quantities of galaxies in a highdimensional space. Various (logarithmic) linear relations, so-called galactic scaling relations, havebeen discovered. Research to unify the scaling relations and find the fundamental relationships has
led to the concept of galaxy manifolds. However, the galaxy manifold has once been almost forgotten
because the classical PCA could treat only linear relations, and it remained a limited concept, though
they are still useful for exploring (log)linear relations of galaxies.
Recently, we discovered a galaxy manifold that expresses the basics of galactic evolution by the
Fisher EM algorithm. Because of its strongly nonlinear spatial structure, it could have never beenfound in previous studies based on the classical PCA. To understand the manifold, a moresophisticated method beyond a mere classification is needed. We focused on a method known as the
manifold learning, one of the latest methods of data science that is completely different fromconventional methodologies
2.2 Galaxy manifold constructed by manifold learning
We adopt the algorithm Isomap and UMAP (Uniform Manifold Approximation and Projection).
Isomap defines the neighboring points by using input-space distance and the distant points as asequence of “short hops” between neighboring points. Isomap tries to find shortest paths in a graphwith edges connecting neighboring data points. By construction, Isomap preserves the “surface
density” of data points in the feature space. UMAP is based on differential geometry and algebraic
topology. The algorithm is founded on three assumptions: 1) the data are uniformly distributed on a
Riemannian manifold, 2) the Riemannian metric is locally constant (or can be approximated as such),and 3) the manifold is locally connected. From these assumptions it is possible to model the manifoldwith a fuzzy topological structure. Sine it defines the manifold so that the data points distribute ashomogeneously as possible, it does not preserve the surface density of data points. UMAP alsopreserves some important structural properties, and it is more robust against noise than Isomap.
Manifold learning algorithm can “unfold” a curved and/or rolled manifold in the feature space, andprovide a local coordinate system on it. The resulting manifolds with local coordinates from Isomap
and UMAP are presented in Fig. 1. From Figure 1, we clearly see that the galaxy manifold is twodimensional. We also stress that two different algorithms, Isomap and UMAP yield similar twodimensional manifolds. The difference of the two estimated manifolds is clearly seen in Fig. 1. Since
Isomap preserves the density of data point cloud, we observe that the manifold has a density structure,i.e., dense and sparse regions on the manifold.
Figure 1: The “unfolded” galaxy manifold by amanifold learning algorithm Isomap and UMAP.
Left and right panels show the manifolds from
Isomap and UMAP, respectively. Though theglobal shape is slightly different from each other,
they share common features on the manifold.

2.3 Result
The galaxy manifold obtained with Isomap preserve this information and reveal the speed of galaxy
evolution at various stages along the manifold. e.g., galaxies passes the green valley very fast. In
contrast, the galaxy manifold obtained with UMAP is imposed uniformity on the galaxy data, leading
to a more robust and representative description of the observed galaxy properties e.g., galaxies evolve
continuously in the feature space, without a discontinuity or “jump” on their evolutionary tracks. Thus,
the galaxy manifold provides a clue to the evolutionary path of galaxies on the manifold. The SFR
and stellar mass fields do not show the same evolutionary path. This supports that the galaxy merger
without star formation plays a significant role in the growth of stellar mass. Next step is to fully
parametrize the evolution equation of galaxies.

Normal-Reference Tests for High-Dimensional Hypothesis Testing
Jin-Ting Zhang
National University of Singapore

Abstract
In the past two decades, much attention has been paid for high-dimensional hypothesis testing.
Several centralized or non-centralized L2-norm based test statistics have been proposed. Most ofthem imposed strong assumptions on the underlying covariance structure of the high-dimensional
data so that the associated test statistics are asymptotically normally distributed. In real data analysis,
however, these assumptions are hardly checked so that the resulting tests have a size control
problem when the required assumptions are not satisfied. To overcome this difficulty, in this talk,
we investigate a so-called normal-reference test which can control the size well. In the normalreference test, the null distribution of a test statistic is approximated with that of a chi-square-type
mixture which is obtained from the test statistic when the null hypothesis holds and when thesamples are normally distributed. The distribution of the chi-square-type mixture can be wellapproximated by a three-cumulant matched χ2-approximation with the approximation parametersconsistently estimated from the data. Two simulation studies demonstrate that in terms of size control,
the proposed normal- reference test performs well regardless of whether

the data are nearly

uncorrelated, moderately correlated, or highly correlated and it performs much better than two
existing competitors. A real data example illustrates the proposed normal-reference test.

KEY WORDS: χ2-type mixtures; high-dimensional data; three-cumulant matched χ2-approximation;
two-sample Behrens–Fisher problem.

Generalized information criterion for
high-dimensional PCA rank selection
Su-Yun Huang(Institute of Statistical Science, Academia Sinica)

Abstract
Principal component analysis (PCA) is a commonly used statistical tool for dimension
reduction. An important issue in PCA is to determine the rank, which is the numberof dominant eigenvalues of the covariance matrix. Among information-based criteria,
the Akaike information criterion (AIC) and the Bayesian information criterion (BIC)
are the two most common ones. Both use the number of free parameters for assessing
model complexity, which requires the validity of the simple spiked covariance model. As aresult, AIC and BIC may suﬀer from the problem of model misspeciﬁcation when the taileigenvalues do not follow the simple spiked model assumption. To alleviate this diﬃculty,
we adopt the idea of the generalized information criterion (GIC) to propose a model
complexity measure for PCA rank selection. The proposed model complexity takes into
account the sizes of eigenvalues and, hence, is more robust to model misspeciﬁcation.
Asymptotic properties of our GIC are established under the high-dimensional setting,where n goes to inﬁnity and p/n goes to a constant c > 0. Our asymptotic results showthat GIC is better than AIC in excluding noise eigenvalues, and is more sensitive than
BIC in detecting signal eigenvalues. Numerical examples will be presented.

1