Matey Neykov

Boris P Hejblum

Jennifer A Sinnott
email: jsinnott@stat.osu.edu

Kernel Machine Score Test for Pathway Analysis in the Presence of Semi-Competing Risks

Keywords: Kernel Machines, Pathway Analysis, Resampling, Score Test, Semi-Competing Risks

In cancer studies, patients often experience two different types of events: a non-terminal event such as recurrence or metastasis, and a terminal event such as cancer-specific death. Identifying pathways and networks of genes associated with one or both of these events is an important step in understanding disease development and targeting new biological processes for potential intervention. These correlated outcomes are commonly dealt with by modeling progression-free survival, where the event time is the minimum between the times of recurrence and death. However, identifying pathways only associated with progression-free survival may miss out on pathways that affect time to recurrence but not death, or vice versa. We propose a combined testing procedure for a pathway's association with both the cause-specific hazard of recurrence and the marginal hazard of death. The dependency between the two outcomes is accounted for through perturbation resampling to approximate the test's null distribution, without any further assumption on the nature of the dependency. Even complex nonlinear relationships between pathways and disease progression or death can be uncovered thanks to a flexible kernel machine framework. The superior statistical power of our approach is demonstrated in numerical studies and in a gene expression study of breast cancer.

Introduction

Identifying genomic disruptions associated with cancer prognosis is an important first step in improving scientific understanding of the progression of the disease and identifying important biological processes to potentially target for treatment development. While we could pursue this goal gene-by-gene, interrogating data at the pathway level may be a more appealing approach. Examining pathways leverages the extensive and ever-expanding research of biologists working to improve our knowledge of how genes work together in groups. For example, the Molecular Signatures Database 1 now contains eight major collections of pathways, including a collection of so-called hallmark gene sets which culls information on fifty important biological states and processes from founder gene sets proposed by many different research groups. Moreover, there has been evidence in certain cancers that while many individuals will experience dysregulation in a few key pathways, the actual component of the pathway that is altered may differ across individuals. 2;3 Thus, focusing on commonly dysregulated pathways instead of commonly dysregulated genes may be a more effective approach for identifying key processes in disease progression.

A number of statistical methods have emerged in recent years to assess the significance of the association between a set of genes and an outcome of interest [START_REF] Goeman | A global test for groups of genes: testing association with a clinical outcome[END_REF][START_REF] Goeman | Testing association of a pathway with survival using gene expression data[END_REF][START_REF] Liu | Semiparametric regression of multidimensional genetic pathway data: Least-squares kernel machines and linear mixed models[END_REF][START_REF] Liu | Estimation and testing for the effect of a genetic pathway on a disease outcome using logistic kernel machine regression via logistic mixed models[END_REF][START_REF] Cai | Kernel Machine Approach to Testing the Significance of Multiple Genetic Markers for Risk Prediction[END_REF][START_REF] Lin | Kernel machine SNP-set analysis for censored survival outcomes in genome-wide association studies[END_REF][START_REF] Wu | Rare-variant association testing for sequencing data with the sequence kernel association test[END_REF][START_REF] Wu | Kernel machine snp-set testing under multiple candidate kernels[END_REF][START_REF] Sinnott | Omnibus risk assessment via accelerated failure time kernel machine modeling[END_REF][START_REF] Hejblum | Time-Course Gene Set Analysis for Longitudinal Gene Expression Data[END_REF] . Methods based on kernel machines (KMs) are particularly attractive because they allow a lot of modeling flexibility while relying on specification of "similarity" between two individuals' pathway expression, which can be easier to conceptualize and discuss with disease researchers. KM methods can be implemented with a linear kernel to induce a model with linear effects, or with a variety of nonlinear kernels to allow increasingly complex models with nonlinear effects without explicitly specifying the form of those effects. Furthermore, when the sampling distributions of KM tests are estimated with perturbation resampling, it is easy to combine information across tests using different kernels and control for multiple testing when many pathways are under consideration. 8;11;12 .

When testing for a pathway's association with survival, methods have been developed for the Cox proportional hazards model and the accelerated failure time model. 8;12 However, in studies of disease survival, it is common for patients to experience both non-terminal and terminal disease related events. For example, in studies of survival following breast cancer diagnosis, patients may experience both recurrence of the disease (a non-terminal event) and death (a terminal event). One common approach for dealing with this kind of data is simply to ignore information on the non-terminal event and focus on overall survival (OS); this approach may be appealing because OS is the most clinically relevant and accurately measured endpoint. 14;15 A second common approach is to model so-called progressionfree survival (PFS), in which the failure time is the time to either recurrence or death -whichever comes first. [START_REF] Broglio | Detecting an overall survival benefit that is derived from progression-free survival[END_REF][START_REF] Buyse | Survival is not a good outcome for randomized trials with effective subsequent therapies[END_REF][START_REF] Amir | Poor correlation between progression-free and overall survival in modern clinical trials: are composite endpoints the answer?[END_REF] . Using PFS has some benefits, especially in clinical trials. In clinical trials, focusing on PFS allows assessment of treatment efficacy sooner, and thus is more efficient when progression is a surrogate marker for death. [START_REF] Sherrill | Review of meta-analyses evaluating surrogate endpoints for overall survival in oncology[END_REF] A drug's impact on PFS may also be a more accurate assessment of its value because patients' treatment plans are likely reassessed at the time of recurrence without regard to randomization. 16;17 However, in a study relating tumor gene expression to survival, methods that use only PFS may lose substantial information by eliminating follow-up time contributed by patients after they experience recurrence. Moreover, it is plausible that a biological process hastens time to recurrence but not time to death, or vice versa; or that dysregulation in some subnetwork of a pathway brings about recurrence while dysregulation in another accelerates the time to death. Identifying such biological processes and pathways may be important for both prediction and biological understanding, but may be difficult if analysis focuses only on a composite endpoint combining recurrence and death. Indeed, our simulations demonstrate that methods based on PFS can have poor performance in exactly these situations.

To more fully capture disease progression integrating all available information on recurrence and death, numerous methods are emerging to address this so-called Semi-Competing Risks (SCR) setting. The main statistical challenge for SCR data is that the marginal survival function for the non-terminal event is not nonparametrically identifiable. [START_REF] Fine | On semi-competing risks data[END_REF] Specifically, let T R be the time to the non-terminal event and T D be the time to the terminal event (R and D here stand for "recurrence" and "death"). Then we can define the joint survival function S(t R , t D ) = P (T R > t R , T D > t D ), but we can only learn about this function on the "wedge" region 0 < t R < t D , because recurrence times after death are never observable. Because of this, the marginal survival function of T R is not estimable without further modeling assumptions. Some methods to deal with SCR model the two marginal distributions of T R and T D , and either leave their joint distribution arbitrary or model it using a copula. [START_REF] Fine | On semi-competing risks data[END_REF][START_REF] Cook | Marginal analysis of recurrent events and a terminating event[END_REF][START_REF] Ghosh | Nonparametric analysis of recurrent events and death[END_REF][START_REF] Ghosh | Marginal regression models for recurrent and terminal events[END_REF][START_REF] Jiang | Semiparametric analysis of survival data with left truncation and dependent right censoring[END_REF][START_REF] Ghosh | Semiparametric inferences for association with semi-competing risks data[END_REF][START_REF] Peng | Regression modeling of semicompeting risks data[END_REF][START_REF] Lakhal | Estimating survival and association in a semicompeting risks model[END_REF][START_REF] Hsieh | Regression analysis based on semicompeting risks data[END_REF][START_REF] Fu | Joint modeling of progression-free survival and overall survival by a bayesian normal induced copula estimation model[END_REF] Other methods model the two cause-specific hazards for recurrence and death, as well as a third conditional hazard for the time to death following recurrence. [START_REF] Liu | Shared frailty models for recurrent events and a terminal event[END_REF][START_REF] Ye | Semiparametric analysis of correlated recurrent and terminal events[END_REF][START_REF] Zeng | Semiparametric transformation models with random effects for joint analysis of recurrent and terminal events[END_REF][START_REF] Xu | Statistical analysis of illness-death processes and semicompeting risks data[END_REF][START_REF] Zeng | Estimating treatment effects with treatment switching via semicompeting risks models: an application to a colorectal cancer study[END_REF][START_REF] Zhang | Bayesian gamma frailty models for survival data with semi-competing risks and treatment switching[END_REF][START_REF] Hu | Joint modeling approach for semicompeting risks data with missing nonterminal event status[END_REF][START_REF] Lee | Bayesian semiparametric analysis of semicompeting risks data: investigating hospital readmission after a pancreatic cancer diagnosis[END_REF] To our knowledge, previous methods for the SCR setting (as well as the competing risks

[CR] setting) have focused on modeling linear effects of the covariates on the log-hazards, but have not been designed to effectively capture nonlinear effects of covariates on the log-hazards.

To test whether a pathway is associated with recurrence and/or death, we propose a novel approach, in which we separately test whether the pathway is associated with the cause-specific hazard for recurrence Prepared using sagej.cls and the marginal hazard for death. This approach is more appealing than the first set of approaches that model the marginal distributions of T R and T D because it only imposes modeling assumptions on observable quantities. It also provides more flexibility than the second set of approaches that model causespecific and conditional hazards because it only imposes structure on two hazards -the cause-specific hazard of recurrence and the overall hazard of death. The two component pathway tests are KM-based tests, and thus allow flexible, data-driven detection of any nonlinear effects that the pathway genes may have on either the non-terminal or terminal events. In order to combine information from the two tests, we add the test statistics together, which is mathematically equivalent to making a working independence assumption about T R and T D . In actuality, the two event time are likely to be correlated; therefore, we estimate the null distribution of this combined test statistic using a perturbation resampling approach that preserves the correlation between the two event times, but does not explicitly specify the nature of the relationship.

The paper proceeds as follows. In section 2 we describe our KM SCR framework: in particular, in subsection 2.3 we develop a score test extending the work of Cai et al. [START_REF] Cai | Kernel Machine Approach to Testing the Significance of Multiple Genetic Markers for Risk Prediction[END_REF] for the Cox KM model; full details of the derivation as well as the asymptotic behavior of the test are detailed in the Appendix. In subsection 2.4 we propose a perturbation resampling procedure to approximate the null distribution of the score test. In section 3, we present numerical studies, in which we compare our proposed approach to pathway tests for three alternative models commonly used in practice -one modeling PFS, one modeling OS alone, and one modeling the two event types as CR. Our proposed approach is illustrated in an application relating breast tumor gene expression to recurrence and death in section 4, and we demonstrate how our SCR model improves over the more common PFS model in real data. Some final comments are made in section 5. The approach is implemented in the R package kernscr available on CRAN at https://cran.rproject.org/web/packages/kernscr.

Methods

Problem Description and Notation

Let T R be the time to the non-terminal event and T D be the time to the terminal event. In our setting, T R is the time to cancer recurrence, and T D is the time to cancer-specific death. There are four possible observation patterns -we may observe: 1. T R followed by T D : patients who experience the non-terminal event followed by the terminal event 2. T R only: patients who experience the non-terminal event but do not experience the terminal event during follow-up Prepared using sagej.cls 3. T D only: patients who experience the terminal event without the non-terminal event occurring 4. neither T R nor T D : patients who do not experience either kind of event during follow-up We let C be the time of censoring for each subject, which is the time of last follow-up or the time of death due to another cause; C is observed for subjects with patterns 2 and 4. The goal of the study is to evaluate whether a set of genes (or SNPs, or other biomarkers) is associated with recurrence and/or death. The genes in the gene set will be denoted by Z p×1 , a vector of length p. To have the best chance of identifying a gene set that is relevant to progression and can add value to existing clinical predictors, we may want to include information on known clinical variables U such as tumor grade and stage. We will assume that the variable U is discrete and is coded to take values {1, 2, . . . , K} -conceptually, these represent different risk strata, and will be used as strata in our model. For example, in breast cancer, the levels of U could correspond to different risk strata defined by combinations of factors in the Gail model. [START_REF] Gail | Projecting individualized probabilities of developing breast cancer for white females who are being examined annually[END_REF] Because T D may censor T R , and C may censor both times, we define

X D = T D ∧ C, ∆ D = I(T D ≤ C), X R = T R ∧ X D , and ∆ R = I(T R ≤ X D ).
Here a ∧ b denotes the minimum of a and b, and

I(A)
is the indicator function taking the value 1 if A occurs and 0 otherwise. For convenience of notation, we also define W = (Z T , U ) T . The observed data for n iid individuals is then:

O = X R i , X D i , ∆ R i , ∆ D i , W T i : i = 1, . . . , n .
We assume that censoring is non-informative and that censoring is independent of the event times conditional on the covariates: C ⊥ (T R , T D ) | W. The time to death T D may also censor T R , but these two times are likely to be correlated, so that this censoring is informative. We will describe several approaches that could be used to model this data. These approaches will impose structure on either the cause-specific hazard of a certain event E ∈ {R, D} -which we will denote by λ c -or the marginal hazard of an event E ∈ {R, D} -which we will denote by λ m . That is, we will describe different approaches that model either the instantaneous rate of occurrence of event E given that the individual is still at risk of experiencing both R and D -the cause-specific hazard:

λ E c (t; W) = lim ε→0 P (t ≤ T E < t + ε | T R ≥ t, T D ≥ t, W)/ε, E ∈ {R, D}
or the instantaneous rate of occurrence of event E given that the individual is still at risk of E, without regard to the other event type -the marginal hazard:

λ E m (t; W) = lim ε→0 P (t ≤ T E < t + ε | T E ≥ t, W)/ε, E ∈ {R, D} Prepared using sagej.cls

The Kernel Machine Semi-Competing Risks Model

Our proposed testing method will be abbreviated "SCR" to indicate Semi-Competing Risks. We model the cause-specific hazard for the non-terminal event and the marginal hazard for the terminal event:

λ R c (t; W) = λ R cU (t) exp {h R (Z)} λ D m (t; W) = λ D mU (t) exp {h D (Z)} (1)
Here, the cause-specific hazard λ R c (t; W) that a person with covariates W experiences recurrence is the product of an unspecified cause-specific baseline hazard λ R cU (t) associated with the particular clinical risk stratum defined by U and the exponential of some function h R of the genomic covariates. The marginal hazard of death for the same person is modeled similarly. The functions h R and h D are unknown, centered, smooth functions of the genomic covariates, and we will ultimately be interested in testing H 0 : h R (•) = 0, h D (•) = 0, within risk strata defined by U. This null suggests that none of the genes in the pathway is associated with either recurrence or death. Note that under the null, the model makes no proportional hazards assumptions, assuming only strata-specific unspecified hazard functions; we discuss this further in the Discussion (Section 5).

We propose to specify the functions h R and h D using the KM framework. Specifically, we assume that both h R and h D belong to H K , the Hilbert space generated by a chosen positive definite kernel K(•, •; ρ). The kernel K is a measure of similarity between two vectors of genomic measurements, and may depend on a possibly unknown scaling parameter ρ. Different choices of kernel K will yield different collections of possible functions h(•). For example, specifying the linear kernel K(z 1 , z 2 ) = z T 1 z 2 forces h(z) = β T z, a linear function of the covariates. The quadratic kernel K(z 1 , z 2 ; ρ) = (ρ + z T 1 z 2 ) 2 yields a Hilbert space H K spanned by basis functions {z j , z j z j : j, j = 1, . . . , p}, which incorporates main effects, quadratic effects, and 2-way interactions. To allow for more complex non-linear effects, one may consider the Gaussian kernel, defined by K(z 1 , z 2 ; ρ) = exp{-z 1 -z 2 2 /ρ}. The resulting function space H K is generated by the radial basis functions. Many other kernels are available for use depending on the data available and the assumptions researchers wish to make: these include the Bessel kernel, the B-Spline kernel, and the sigmoid kernel. Adaptations such as kernel principal components analysis (PCA) can be implemented to reduce the complexity of the feature space associated with the kernel function and thus further improve the performance of KM methods in data sets of moderate sample size. 8;39 When we do not assume the effects of the pathway variables Z on recurrence and death are identically 0, we may obtain estimators for each h E (for E ∈ {R, D}) by maximizing the associated partial Prepared using sagej.cls likelihood functions each penalized by the norm of h E in the Hilbert space, h E H K :

log(L E P (h E )) = n i=1 K k=1 I(U i = k) log exp{h E (Z i )} n l=1 I(U l = k)Y E l (s) exp{h E (Z l )} dN E i (s) - d E 2 h E 2 H K .
where

N E i (t) = ∆ E i I(X E i ≤ t) for both event types, but the at-risk indicators differ: Y R i (t) = I(X R i ≥ t)I(X D i ≥ t) and Y D i (t) = I(X D i ≥ t).
The d E are penalty parameters that control the smoothness of the functions h E . By the representer theorem, 40;41 the maximizer of log(

L E P (h E )) for h E must take a dual form, h E (z) = n l=1 α El K(Z l , z; ρ)
, where the α El are unknown parameters. Under this dual representation, we maximize (with a slight abuse of notation):

log(L E P (α E )) = n i=1 K k=1 I(U i = k) log exp(α E K i (ρ)) n l=1 I(U l = k)Y E l (s) exp(α E K l (ρ)) dN E i (s) - d E 2 α E K(ρ)α E (2)
where

K i (ρ) = (K i1 (ρ), K i2 (ρ), . . . K in (ρ))
is the i th row of the kernel matrix K(ρ). The penalized partial likelihood functions (2) may be viewed as arising from the random effects model in which we assume:

λ R ci (t) = λ R cUi (t) exp {α T R K i (ρ)} λ D mi (t) = λ D mUi (t) exp {α T D K i (ρ)} (3)
for

α E = τ E ε E , where ε E is multivariate normal with E[ε E ] = 0 and Var(ε E ) = K(ρ) -, and d -1 E = τ 2 E .
Here K(ρ) -is the Moore-Penrose generalized inverse of K(ρ). Analogous connections between KM models and the mixed model framework have successfully been used to fit other KM regression models. 6-8;12 Conditional on the random effects ε E , the partial likelihoods may be written:

log(L E ε E (τ E ; ρ)) = n i=1 K k=1 I(U i = k) log exp{τ E ε T E K i (ρ)} U (0) Ek (s; ρ, τ E ) dN E i (s) (4)
where we define

U (m) Ek (s; ρ, τ E ) = n i=1 I(U i = k)Y E i (s) exp{τ E ε T E K i (ρ)} • {ε T E K i (ρ)} m for m = 0, 1, 2.

Score test

With the model formulated as in (3), the null hypothesis becomes H 0 : τ R = 0, τ D = 0. Because these values are on the boundary of the parameter space, the likelihood ratio test has a nonstandard distribution; however, we can derive a score test similar to that derived in Cai et al. [START_REF] Cai | Kernel Machine Approach to Testing the Significance of Multiple Genetic Markers for Risk Prediction[END_REF] for the Cox KM model with a Prepared using sagej.cls single event type of interest. To proceed in our setting, we will derive a test statistic under an assumption of working independence for T R and T D . This gives us a guide for how to combine information about the relationship between the covariates and each event time, but we would not expect these two times to be independent in practice. Thus, although the working independence assumption is used to justify the choice of the test statistic, both the asymptotic distribution we derive in the Appendix (Section 5) and the resampling approach used to estimate the null distribution in finite sample (Subsection 2.4) are valid when the two times are correlated, as we expect to occur in practice.

If we assume working independence for T R and T D ; independence of the random effects ε R and ε D ; and that we are testing H 0 against an alternative of the form H A : τ R = τ D η = τ for some known constant η; then the joint likelihood is:

M(τ ; ρ) = log L R ε R (τ ; ρ) + log L D ε D (ητ ; ρ) .
Then, following Cai et al. [START_REF] Cai | Kernel Machine Approach to Testing the Significance of Multiple Genetic Markers for Risk Prediction[END_REF] and Commenges and Andersen, [START_REF] Commenges | Score Test of Homogeneity for Survival Data[END_REF] a score test may be based on the statistic:

Q(ρ, η) = E ε R ,ε D ∂M(τ ; ρ) ∂τ τ =0 2 + E ε R ,ε D ∂ 2 M(τ ; ρ) ∂τ 2 τ =0
.

Derivations are in the Appendix. In the case where η = 1 the form of the score statistic

Q(ρ) = Q(ρ, 1)
is:

Q(ρ) = E∈{R,D} Q E (ρ) where Q E (ρ) = M E (∞) K(ρ) M E (∞) -n q E (ρ) (5)
where

M E (s) = ( M E 1 (s), M E 2 (s), . . . , M E n (s)) T and each M E i (s) = N E i (s) - s 0 Y E i (u)d Λ E Ui (u), Λ E k (s) = n j=1 I(U j = k) s 0 dN E j (u) S (0) Ek (u) , S (0)
Ek (s) = n i=1 I(U i = k)Y E i (s), and q E (ρ) = 1 n K k=1 n j=1 I(U j = k)Y E j (s)K jj (ρ)d Λ E k (s) - 1 n K k=1 n j=1 n l=1 I(Uj =k)I(U l =k)Y E j (s)Y E l (s)K jl (ρ) S (0) Ek (s) d Λ E k (s).
Prepared using sagej.cls

A resampling procedure for approximating the null distribution

In the Appendix, we show that n -1 Q(ρ) + q(t) is asymptotically equivalent to the process:

W(ρ) = E∈{R,D} K(z 1 , z 2 ; ρ)d W E M (z 1 )d W E M (z 2 ) + -2 K k=1 p k ω E k (z, t, ρ)d W E Λ E k (t) d W E M (z) + K k=1 K k =1 p k p k ω E kk (t, s, ρ)d W E Λ E k (t)d W E Λ E k (s)
where

W E M (z) = n -1 2 n i=1 M E i (∞)I(Z i ≤ z), M E i (t) = N E i (t) - t 0 Y E i (s)dΛ E Ui (s), p k = P (U = k), ω E k (z, t, ρ) = E[K(Z, z; ρ)Y E (t) | U = k], ω E kk (t, s, ρ) = E[K(Z 1 , Z 2 ; ρ)Y E 1 (t)Y E 2 (s) | U 1 = k, U 2 = k ].
This form is very close to the form given in Cai et al.; 8 the main differences are that we are summing the test statistics over the two event types, and that the covariates U are modeled as strata.

Also note that Λ R k (t) = Λ R ck (t) and Λ D k (t) = Λ D mk (t)
, but we drop the subscripts distinguishing the types of hazards being modeled for notational simplicity. As discussed in the Appendix, the arguments in Cai et al. [START_REF] Cai | Kernel Machine Approach to Testing the Significance of Multiple Genetic Markers for Risk Prediction[END_REF] may be easily generalized to our setting, in order to establish that W(ρ) converges weakly to a zero-mean Gaussian process W(ρ).

In finite samples, we may approximate the distribution of Q(ρ) by perturbation resampling.

Specifically, suppose we generate a vector of iid standard normal variates V = (V 1 , . . . , V n ) T , and calculate:

W * (ρ) = E∈{R,D} K(z 1 , z 2 ; ρ)d W E * M (z 1 )d W E * M (z 2 ) + -2 K k=1 p k ω E k (z, t, ρ)d W E * Λ E k d W E * M (z) + K k=1 K k =1 p k p k ω E kk (t, s, ρ)d W E * Λ E k (t)d W E * Λ E k (s)
where

W E * M (z) = n -1 2 n i=1 M E i (∞)I(Z i ≤ z)V i , W E * Λ E k (t) = n -1 2 n i=1 W E ki (t)V i , ω E k (z, t, ρ) = 1 n n i=1 I(U i = k)K(Z i , z; ρ)Y E i (t), ω E kk (t, s, ρ) = 1 n 2 n i=1 n j=1 I(U i = k)I(U j = k )K(z i , z j ; ρ)Y E i (t)Y E j (s), p k = 1 n n i=1 I(U i = k), and where W E ki (t) = t 0 I(Ui=k)d M E i (s) S (0) Ek (s)
is the empirical counterpart of W E ki (t), defined in the Appendix. By assigning each individual a perturbation variate V i and associating this variate with the individual in each process

Prepared using sagej.cls W R * M (z), W R * Λ R k (t), W D * M (z), W D * Λ R k (t),
-1 Q(ρ) by Q * (b) (ρ) = W * (b) (ρ) -q * (ρ)
, where q * (ρ) is the empirical mean of W * (ρ). For example, for any fixed ρ, an approximate p-value for the test may be calculated by 1

B B b=1 I Q * (b) (ρ) > n -1 Q(ρ) .

Kernels with tuning parameters

The linear kernel does not rely on a tuning parameter ρ, but many kernels, including the quadratic and Gaussian kernels used here, rely on a tuning parameter that controls the complexity of the associated functions h. Under H 0 , the matrix K(ρ) disappears, so that ρ is not estimable; however, we may follow the approach in Cai et al. [START_REF] Cai | Kernel Machine Approach to Testing the Significance of Multiple Genetic Markers for Risk Prediction[END_REF] of using the statistic:

S = sup ρ∈[ρ L ,ρ U ] n -1 Q(ρ)/ σ(ρ)
where σ 2 (ρ) is an estimate of the variance σ 2 (ρ) of W(ρ). To get such an estimate, we can calculate the variance of the perturbations: σ 2 (ρ) = Var W * (ρ) . We then may calculate an approximate p-value by calculating 1 For the Gaussian and quadratic kernels, we determine the range [ρ L , ρ U ] of ρ by requiring that the associated kernel matrices K(ρ) have eigenvalues λ i that decay at a polynomial rate O(i -α ) for some range of α > 1. This approach is motivated by work by Braun, [START_REF] Braun | Spectral properties of the kernel matrix and their application to kernel methods in machine learning[END_REF] who bounds the error due to projecting the feature space onto the first r principal components using terms whose behavior depends on the decay properties of the eigenvalues. By experimentation we found that α ∈ [1.5, 4] yielded a reasonable range of feature space complexity. To implement this, for each ρ, we regressed the logarithms of the top 90% of eigenvalues on their indices and used the regression coefficient as an estimate of α; we considered values ρ whose associated α fell in the specified range.

Alternative Models

As discussed in the introduction, a popular modeling alternative for this setting is to focus on progression-free survival (PFS). To do this, let RD stand for recurrence-or-death, and define

T RD = T R ∧ T D , X RD = T RD ∧ C, ∆ RD = I(T RD ≤ C).
Then we model the (necessarily marginal) hazard Prepared using sagej.cls of recurrence-or-death by:

λ RD m (t) = λ RD mU (t) exp {h RD (Z)} ,
where h RD (•) ∈ H K for some kernel K. The previous method developed in Cai et al. [START_REF] Cai | Kernel Machine Approach to Testing the Significance of Multiple Genetic Markers for Risk Prediction[END_REF] applies directly.

Another popular modeling approach is to disregard information on metastases, and focus entirely on the marginal model for death -we will call this model OS for overall survival:

λ D m (t; W) = λ D mU (t) exp {h D (Z)} .
Once again, the previous method developed in Cai et al. [START_REF] Cai | Kernel Machine Approach to Testing the Significance of Multiple Genetic Markers for Risk Prediction[END_REF] applies directly.

Another alternative model which we will compare is the competing risks model (CR), in which we model both cause-specific hazards of the times to recurrence and death and disregard follow-up time after recurrence occurs:

λ R c (t; W) = λ R cU (t) exp {h R (Z)} λ D c (t; W) = λ D cU (t) exp {h D (Z)}
A score test for this model may be derived in the same manner as outlined above for the SCR model, under an assumption of working independence for the two event times. By comparing the power of this model to the power of the SCR model, we can quantify the added benefit seen by not just separating recurrence and death, but also by including information on death after recurrence.

Simulation study

We present the results of simulation studies comparing KM-based testing under our proposed SCR model and the three alternative models described in Subsection 2.6 -the PFS model, the CR model, and the OS model. We used three kernels -linear, quadratic, and Gaussian. All tests were run using 1000 perturbations with perturbation variables V generated from a standard normal distribution. The numbers presented are based on 2000 simulations for empirical size and 1000 simulations for empirical power.

For simplicity, we do not consider any risk stratifying clinical covariate U . For p = 5 and 10, we generate genomic covariates Z p×1 independently for each individual, where Z is multivariate normal with mean 0, variance 1, and correlation between each pair of genes given by ρ Z . We consider ρ Z = 0.2 to reflect low correlation between genes and ρ Z = 0.5 to reflect moderate correlation between genes. We consider studies of size n = 100 and n = 200. The times until the non-terminal event and the terminal Prepared using sagej.cls event are generated according to:

T R = exp{C R + h R (Z) + E R } T D = exp{C D + h D (Z) + E D }.
Here, the error terms E R and E D marginally follow extreme value distributions but are generated to be correlated. Specifically, we let Y ∼ MVN((0, 0) T , Σ), where Σ = ( 1 0.5 0.5 1 ); then we transform E = log{-log{Φ(Y)}}, where Φ is the standard normal cdf, and let the first value be E R and the second value be E D . The functions h R (Z) and h D (Z) differ across simulation settings, and the constants C R and C D are selected to adjust the relative sizes of T R , T D , and C, which is the censoring time generated independently from an exponential distribution with mean 1 30 . Thus, correlation between the recurrence time T R and the death time T D is induced both through the error terms E R and E D , and through h R and h D if they are nonzero and involve common or correlated genes.

Table 1 displays the empirical sizes across different data simulation settings, for each of the four models and each of the three kernels. In this setting, approximately 31% recur and then die; 14% only recur during follow-up; 28% die (without recurrence) before the end of follow-up, and 27% are censored before they experience either type of event. In general, the sizes appear to be well-calibrated. The tests using the linear kernel have size near 5%, while the tests using the nonlinear kernels have levels below that. The tests using the Gaussian and quadratic kernels are particularly conservative when the betweengene correlation is low (ρ Z = 0.2) and n = 100. On average, the test levels improve with increasing sample size n.

In Table 2, we compare test performance when the true signal is linear. Specifically, we let the signal for recurrence be h R (Z) = Z 1 /4, and the signal for death be h D (Z) = Z 2 /4 -Z 3 /4 + Z 4 /4. In this setting, approximately 30% experience both recurrence and death; 14% experience recurrence but not death during follow-up; 28% experience death without recurrence; and 28% are censored. Here, not surprisingly, the linear kernel (L) typically performs the best for all the models. In this setting, we expect the SCR model to outperform the PFS model because the PFS model merges the two event types, but the genes that are associated with time to death and time to recurrence differ (though they are correlated).

And indeed, under the linear kernel, the SCR model typically outperforms the others -for example, when n = 200, p = 5 and ρ Z = 0.5, and the linear kernel is used, the power to detect the pathway with the SCR-based test is 77.1%. In comparison, the PFS-based test has power of only 62.7%. The CR-based test has an even lower power of 60.8%, thus demonstrating a gain in power of 16.3% by the SCR test from including information after recurrence for the 44% of subjects who recur during follow-up. Finally, the OS test has power 67.8%, reflecting some decrease in power when information on the recurrence time is not taken into account -though interestingly, in this setting the OS-based test is frequently a strong Prepared using sagej.cls Table 1. Empirical sizes (in %) -i.e., the proportion of significant p-values when both hR(•) ≡ 0 and hD(•) ≡ 0 -of the four tests under consideration, at level α = 5%. The tests compared are: SCR (Semi-competing Risks); PFS (Progression-Free Survival); CR (Competing Risks); and OS (Overall Survival). For each of these model formulations, KM (kernel machine) tests are run with Gaussian (G), Linear (L), and Quadratic (Q) kernels. Comparisons are made for sample sizes n = 100 and n = 200; number of pathway variables p = 5 and p = 10; and pathway variable pairwise correlation ρ Z = 0.2 and 0.5. In this setting, approximately 31% recur and then die; 14% only recur during follow-up; 28% die (without recurrence) before the end of follow-up, and 27% are censored before they experience either type of event. performer, and is sometimes more powerful when nonlinear kernels are used, perhaps due to the strength of the signal h D (Z) and the involvement of multiple genes.

n = 100 n = 200 SCR PFS CR OS SCR PFS CR OS p = 5 ρ = 0.2 G 2.2 2.
We also consider a nonlinear setting with h R (Z) = 1.5 sin(Z 1 + Z 2 2 ), h D (Z) = (Z 4 -Z 5 ) 2 /4 (Table 3). Here, approximately 53% of subjects experience both recurrence and death; 20% experience recurrence without death; 18% die before recurring; and 9% are censored. Once more the SCR model typically outperforms the others, and the Gaussian kernel typically outperforms the linear and even the quadratic kernels. The OS model is particularly weak in this setting because it ignores the additional information provided by the 73% of subjects who experience recurrence. For example, when n = 200, p = 5, and ρ Z = 0.5, and the Gaussian kernel is used, the power to detect the pathway is 95.1% with the SCR model; 88.2% when the PFS model is used; 94.5% when the CR model is used; and only 9.2% when OS is used. When n = 100, p = 10, and ρ Z = 0.2, the Gaussian kernel has slightly lower power than the linear kernel. This is likely because of the difficulty in finding the best tuning parameter for the nonlinear kernel when there is not much data, the covariates are only weakly correlated, and the signal in the pathway is sparser. We see this phenomenon diminish when we either increase the sample size or increase the correlation among the genomic covariates.

Application to Breast Cancer Gene Expression

Genomics has proven to be of great importance in breast cancer. Besides high risk mutations identified on BRCA1 and BRCA2 genes, numerous genomic biomarkers have been developed for clinical use in this disease. Such discoveries have translated into better prevention, better diagnoses and better treatment for patients. However, further understanding of the molecular pathways of pathogenesis in breast cancer remains critical for the development of better targeted therapies. [START_REF] Olopade | Advances in breast cancer: Pathways to personalized medicine[END_REF] Here, following Cai et al. [START_REF] Cai | Kernel Machine Approach to Testing the Significance of Multiple Genetic Markers for Risk Prediction[END_REF] , we focus on 70 genetic pathways (with numbers of genes ranging from 2 to 235 [median=22]) from the Molecular Signatures Database 1 and their potential association with breast cancer survival.

Cai et al. identified pathways associated with OS regardless of the development of metastases. Here, we compare using the PFS model, the most commonly chosen model for integrating the two types of events, to the proposed SCR model. We apply both the SCR and PFS models using a Gaussian kernel, with kernel PCA [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF] including principal components accounting for 90% of the total variability, to assess the overall effect of each of these pathways on breast cancer progression and survival, using data from a study by van de Vijver et al. [START_REF] Van De Vijver | A gene-expression signature as a predictor of survival in breast cancer[END_REF] Gene expression was measured for 260 patients with primary breast carcinomas from the Netherlands Cancer Institute. The median follow-up time was 8.8 years; 23% died Prepared using sagej.cls Table 3. Power (in %) of detecting a nonlinear signal, specifically letting hR(Z) = 1.5 sin(Z1 + Z 2

2 ) for recurrence and hD(Z) = (Z4 -Z5) 2 /4 for death. The tests are run at α = 5%, and the four compared are as before: SCR (Semi-competing Risks); PFS (Progression-Free Survival); CR (Competing Risks); and OS (Overall Survival). For each of these model formulations, KM (kernel machine) tests are run with Gaussian (G), Linear (L), and Quadratic (Q) kernels. Comparisons are made for sample sizes n = 100 and n = 200; number of pathway variables p = 5 and p = 10; and pathway variable pairwise correlation ρ Z = 0.2 and 0.5. In this setting, approximately 53% recur and then die; 20% only recur during follow-up; 18% die (without recurrence) before the end of follow-up, and 9% are censored before they experience either type of event. after presenting metastases, 9% were censored with metastases, 2% died without metastases, and 66%

were censored before experiencing either event. Because we are testing 70 pathways simultaneously, we use the multiple testing adjustment method using perturbation described in Cai et al. [START_REF] Cai | Kernel Machine Approach to Testing the Significance of Multiple Genetic Markers for Risk Prediction[END_REF] In the SCR model, incidence of metastases was the non-terminal event, and death was the terminal event. PFS was defined as time to either death or metastases, whichever occurred first. Tests were run with 10,000 perturbations. Note that the number of perturbations determines that smallest size p-value that one can reliably estimate; with 10,000 perturbations, the smallest nonzero nominal p-value is 0.0001, and any more significant associations report a p-value of 0. The smallest nonzero adjusted p-value one can reliably estimate is approximately 0.0001 × 70 = 0.007. Smaller p-values could be estimated by running more perturbations. family-wise error rate of 0.05, 24 pathways are significant using the SCR model while 21 are significant using PFS. In this data set, the difference between these two methods is not dramatic, perhaps because

-3 -2 -1 0 < -4 log10 adjusted p-value -1 0 -2 < -2.2
Figure 1. log10 p-values for testing the overall effect of 70 genetic pathways on breast cancer survival based on the KM (kernel machine) score test with Gaussian kernel using kernel PCA (principal components analysis; with a 90% of total variability threshold). The crosses and the squares represent PFS and SCR models, respectively. These results are based on B = 20,000 perturbations.

Prepared using sagej.cls the total proportion of either type of event is not very large. We would expect the PFS-and SCR-based approaches to produce similar results when, for example, a few key genes in a pathway hasten both time to progression and time to death with similar effects. The SCR-based model will be better than the PFS-based model at detecting pathways in which progression and death are associated with different genes' dysregulation within the pathways. And indeed, in this example, the SCR model identifies three more pathways than the PFS model at the family-wise error rate of 0.05. These pathways include the TGFβ pathway and the MAPK pathway, which have been implicated in breast cancer. For example, hyperactivation of the MAPK pathway is directly associated with the appearance of certain subtypes of breast cancer. [START_REF] Oh | Hyperactivation of MAPK induces loss of ERalpha expression in breast cancer cells[END_REF] Similarly, the TGFβ pathway has recently been identified as a good prognostic indicator for certain subtypes of breast cancer, [START_REF] Teschendorff | Improved prognostic classification of breast cancer defined by antagonistic activation patterns of immune response pathway modules[END_REF] while its implication in breast cancer has been known and studied for more than two decades. 48;49 Thus, the increased power of the SCR KM model over the PFS KM model may directly translate into identifying more biologically meaningful associations.

Discussion

In this paper we present a flexible test procedure for assessing whether a pathway is associated with recurrence and/or death. The proposed approach is appealing because it tests marginally for the relationship between the pathway and the most clinically relevant and precisely measured outcome (disease-specific death), and augments the test by testing for a potential relationship between the pathway and the non-terminal event. Using information on both events is preferable even when overall survival is the most clinically relevant endpoint, because integrating information on recurrence improves our ability to identify important pathways if the recurrence endpoint is a meaningful landmark between diagnosis and death, and is especially helpful if the survival time after recurrence is variable and governed by patient response to a variety of post-recurrence treatments. [START_REF] Broglio | Detecting an overall survival benefit that is derived from progression-free survival[END_REF] The proposed approach is implemented in the R package kernscr available on CRAN at https://cran.r-project.org/web/packages/kernscr. As mentioned, one reason that PFS is sometimes preferred, especially in the clinical trial setting, is that it is not affected by changes in treatment that occur after a patient presents with a recurrence. A method such as ours that combines information from recurrence and eventual death may elide over this crucial source of heterogeneity. Although we cannot fully investigate this issue in the breast cancer data because we lack information about treatment, our proposed method is still informative in this setting for identifying potentially important pathways -we just need to be careful about the interpretation. Specifically, a pathway may be associated with death either independently of treatment or in a manner mediated by treatment, if the gene expression impacts either treatment selection or treatment response. Such a pathway is likely worth investigating further regardless. Furthermore, other work with nonlinear KM regression has demonstrated that nonlinear kernels can successfully detect signals that are based Prepared using sagej.cls on mixture models -for example, latent classes with differential covariate effects. [START_REF] Sinnott | Omnibus risk assessment via accelerated failure time kernel machine modeling[END_REF] Thus, nonlinear KM methods may in fact be well suited to detecting pathways exhibiting heterogeneous patient effects -though more research is needed to investigate the impacts of different sources of patient heterogeneity on the efficacy of KM regression methods. Some previous work in the SCR setting has focused on modeling marginal associations with the two event types and modeling their joint relationship using a copula model; by comparison, our approach only makes modeling assumptions on observable quantities, and the relationship between the non-terminal and terminal event is left unspecified. Extending KM methods to copula models warrants further research.

Other previous work models cause-specific hazards for the two events as well as a conditional hazard for the terminal event given that the non-terminal event has happened; by comparison, our approach makes fewer modeling assumptions by modeling only the cause-specific hazard for the non-terminal event and the marginal hazard for the terminal event. It would, however, be interesting to extend KM methods to this model as well, particularly if further information, such as treatment changes made at time of recurrence, were available and could be integrated into the model. Previous work in the SCR setting has not been able to capture covariates' nonlinear effects on the two events. By using a flexible KM framework, our approach is able to capture nonlinear effects through the specification of a nonlinear kernel.

A potential issue with our approach is that it is unclear a priori how to combine the tests for the two events because the events are likely to be correlated and the correlation structure is unknown. To resolve this issue, we add the two test statistics together, and model their joint null distribution using a perturbation resampling method that accounts for their correlation without explicitly specifying it.

Adding together the test statistics may not be the optimal way of combining information for the two event types, and future work could consider other possible approaches, such as taking the minimum pvalue of the two tests. 11;12 Further extensions and adaptations of this model are possible. For example, rather than choosing a particular kernel, we could test several candidate kernels (such as a linear and a nonlinear kernel) and combine the information in an omnibus test using the perturbation resampling structure [START_REF] Sinnott | Omnibus risk assessment via accelerated failure time kernel machine modeling[END_REF] -this could be useful in practical settings if researchers are unsure of what similarity metric may be most appropriate for their data and what type of model complexity they wish to consider.

As outlined in the text, our models currently assume that the non-genotypic covariates are discrete and identify "risk strata" which we use as stratifying variables in the proportional hazards models. We made this choice in order to better guarantee test validity in general situations; specifically, we wanted to avoid having the null distribution of the test statistic rely on the assumption that the cause-specific hazard of recurrence and the marginal hazard of death are both proportional in the clinical covariates. We feel that modeling clinical covariates as strata is sufficient in many applications, where combinations of covariates classify individuals into certain risk strata. However, if it is of interest to include non-genotypic covariates Prepared using sagej.cls as linear effects, and the two proportional hazards models specified under the null seem reasonable, then the proposed method is easily extendible to this case, with those covariates included in a similar way as in Cai et al. (2011) [START_REF] Cai | Kernel Machine Approach to Testing the Significance of Multiple Genetic Markers for Risk Prediction[END_REF] .

Our proposed joint test exhibits systematically higher power when compared with other testing approaches when different components of the pathway are associated with recurrence and death, as demonstrated by our simulation studies. In an analysis of a real gene expression study of breast cancer, it also identifies more pathways than the analogous test based only on PFS. The SCR KM method proposed here is a powerful tool that can allow disease researchers to better leverage the rich information in their data and to better identify biologically significant pathways to investigate more in depth in targeted studies.

Supplemental material

Derivation of the Score Test

Using the notation in section 2.3, we have our objective function M(τ

; ρ) = log(L R ε R (τ ; ρ)) + log(L D ε D (ητ ; ρ)) and we need to calculate ∂M(τ ;ρ) ∂τ τ =0 and ∂ 2 M(τ ;ρ) ∂τ 2 τ =0 . First, for E ∈ {R, D} ∂ ∂τ log L E ε E (τ ; ρ) τ =0 = n i=1 K k=1 I(Ui = k) Ki(ρ) T εE - S (1)
Ek (s; ρ) S

Ek (s)

dN E i (s) = n i=1 K k=1 I(Ui = k) dN E i (s) -Y E i (s)d Λ E k (s) Ki(ρ) T εE = n i=1 M E i (∞)Ki(ρ) T εE Prepared using sagej.cls where S (m) Ek (s; ρ) = U (m) Ek (s; ρ, 0) = n i=1 I(Ui = k)Y E i (s){Ki(ρ) T εE} m -though we may write S (m) Ek (s; ρ) = S (m) Ek (s) since it does not depend on ρ; Λ E k (s) = n j=1 I(Uj = k) s 0 dN E j (u) S (0) Ek (u)
; and

M E i (s) = K k=1 I(Ui = k) s 0 dN E i (u) -Y E i (u)d Λ E k (u) = N E i (s) - s 0 Y E i (u)d Λ E U i (u).
Next, we may calculate:

Eε E ∂ log{L E ε E (τ ; ρ)} ∂τ τ =0 2 = n i=1 n j=1 M E i (∞) M E j (∞)Kij(ρ) = M E (∞) T K(ρ) M E (∞)
since Eε E [εEε T E ] = K(ρ) -and Ki(ρ) T K(ρ) -Kj(ρ) = Kij(ρ). Thus, because we assume εR and εD are independent and both have mean 0,

Eε R ,ε D ∂M(τ ; ρ) ∂τ τ =0 2 = M R (∞) T K(ρ) M R (∞) + η 2 M D (∞) T K(ρ) M D (∞) Next, for E ∈ {R, D} ∂ 2 ∂τ 2 L E ε E (τ, ρ) τ =0 = n i=1 K k=1 I(Ui = k) S (1)
Ek (s; ρ) 2 S (0)

Ek (s) 2 - S (2)
Ek (s; ρ) S = -n qR(ρ) -nη 2 qD(ρ)

Then the final expression for the score test statistic is:

Q(ρ; η) = M R (∞) T K(ρ) M R (∞) -n qR(ρ) + η 2 M D (∞) T K(ρ) M D (∞) -n qD(ρ) .

Asymptotic Distribution of the Test Statistic

Our asymptotic results are outlined below, but follow from those presented in Cai et al. [START_REF] Cai | Kernel Machine Approach to Testing the Significance of Multiple Genetic Markers for Risk Prediction[END_REF] To write out the details more explicitly, we assume that the kernel tuning parameter ρ is in some interval [ρL, ρU ] and that the genomic variables Z have bounded support. We assume that XR and XD have a finite support [0, τ ], that the density of C given covariates Z is continuous and bounded, and that the marginal density of C is bounded away from 0 on [0, τ ]. We assume that the kernel K is continuously differentiable with respect to all its arguments.

To find the asymptotic distribution of n -1 Q(ρ), we first note by Fleming and Harrington [START_REF] Fleming | Counting Processes and Survival Analysis[END_REF] ,

sup t∈[0,τ ] Λ R k (t) -Λ R k (t) + Λ D k (t) -Λ D k (t) = oP (1)
Prepared using sagej.cls and that

W E Λ E k (t) = n 1 2 Λ E k (t) -Λ E k (t) = n -1 2 n i=1 W E ki (t) + OP (n -1 2 ),
where W E ki (t) = t 0

I(U i =k)dM E i (s) G Ek (s)
, for G Ek (s) = E Y E (t) | U = k . We may see that

M E i (s) = M E i (s) - n -1 2 Y E i (t)d W E Λ E U i
(t) + Op(n -1 ), so that by mimicking the arguments in Cai et al. [START_REF] Cai | Kernel Machine Approach to Testing the Significance of Multiple Genetic Markers for Risk Prediction[END_REF] , n -1 Q(ρ) + q(ρ) is asymptotically equivalent to the process:

W(ρ) = E∈{R,D} K(z1, z2; ρ)d W E M (z1)d W E M (z2) + -2 K k=1 p k ω E k (z, t, ρ)d W E Λ E k (t) d W E M (z) + K k=1 K k =1 p k p k ω E kk (t, s, ρ)d W E Λ E k (t)d W E Λ E k (s)
which is a zero-mean stochastic process, where:

W E M (z) =n -1 2 n i=1 M E i (∞)I(Zi ≤ z) M E i (t)

  this resampling procedure preserves the correlation between the statistics associated with R and D. Moreover, conditional on the data, W * (ρ) has the same limiting distribution as W(ρ). Thus, by obtaining B realizations of V, and for each of these realizations b = 1, . . . , B calculating a perturbed statistic W * (b) (ρ), we can approximate the null distribution of W(ρ) by W * (ρ) = { W * (b) (ρ), b = 1, . . . , B}. Then we may approximate the null distribution of n

  S * (b) > S)/B, where S * (b) = sup ρ∈[ρ L ,ρ U ] Q * (b) (ρ)/ σ(ρ) , b = 1, . . . , B.

Table 2 .

 2 Power (in %) of detecting a linear signal, specifically letting hR(Z) = Z1/4 for recurrence and hD(Z) = Z2/4 -Z3/4 + Z4/4 for death. The tests are run at α = 5%, and the four compared are as before: SCR (Semi-competing Risks); PFS (Progression-Free Survival); CR (Competing Risks); and OS (Overall Survival). For each of these model formulations, KM (kernel machine) tests are run with Gaussian (G), Linear (L), and Quadratic (Q) kernels. Comparisons are made for sample sizes n = 100 and n = 200; number of pathway variables p = 5 and p = 10; and pathway variable pairwise correlation ρ Z = 0.2 and 0.5. In this setting, approximately 30% recur and then die; 14% only recur during follow-up; 28% die (without recurrence) before the end of follow-up, and 28% are censored before they experience either type of event.

Figure 1

 1 Figure 1 displays the raw p-values and the adjusted p-values for the proposed score test from both models. As expected, the SCR model shows greater power. Indeed p-values from the SCR model are typically lower than those from the PFS model (51 out of 70 unadjusted; 58 out of 70 adjusted). With

2 Sthen

 2 Uj = k)Y E j (s) Kj(ρ) T εE Uj = k)I(U l = k)Y E j (s)Y E l (s) Kj(ρ) T εEε T E K l (ρ) S = k)Y E j (s) Kj(ρ) T εEε T E Kj(ρ) d Λ E k Uj = k)I(U l = k)Y E j (s)Y E l (s)K jl (ρ) S Eε R ,ε D ∂ 2 M(τ ; ρ) ∂τ 2 τ =0

  . Our test statistic is the sum of two statistics of the form presented in Cai et al.[START_REF] Cai | Kernel Machine Approach to Testing the Significance of Multiple Genetic Markers for Risk Prediction[END_REF] ; one arising from T R and one arising from T D . The mainpoint to jusitfy is that W R M (z), W D M (z), W R This holds because the functions I(Z ≤ z), M R i (t), and M D i (t) have finite pseudodimensions. Because of this joint convergence, the asymptotic results from Cai et al.[START_REF] Cai | Kernel Machine Approach to Testing the Significance of Multiple Genetic Markers for Risk Prediction[END_REF] generalize to our setting.

	Λ E 1	(t), . . . , W R Λ E K	(t), W D Λ E 1	(t), . . . , W D K Λ E	(t) converge jointly to zero-
	mean Gaussian processes W R M (z), W D M (z), W R Λ E 1	(t), . . . , W R Λ E K	(t), W D Λ E 1	(t), . . . , W D Λ E

K (t)

.

Prepared using sagej.cls

Acknowledgements BPH was supported by the National Institutes of Health (NIH) grant U54 HG007963. JAS was supported by the National Institutes of Health (NIH) grant T32 CA09001 and the A. David Mazzone Career Development Award.

The authors express their deepest gratitude to Professor Tianxi Cai, Harvard University, for her help in preparing this manuscript.


Nyoman Juniarta
email: nyoman.juniarta@loria.fr

Victor Codocedo
email: victor.codocedo@inria.cl

Miguel Couceiro
email: miguel.couceiro@loria.fr

Amedeo Napoli
email: amedeo.napoli@loria.fr

Biclustering Based on FCA and Partition Pattern Structures for Recommendation Systems

Keywords: biclustering, FCA, pattern structures, recommendation

This paper focuses on item recommendation for visitors in a museum within the framework of European Project CrossCult about cultural heritage. We present a theoretical research work about recommendation using biclustering. Our approach is based on biclustering using FCA and partition pattern structures. First, we recall a previous method of recommendation based on constant-column biclusters. Then, we propose an alternative approach that incorporates an order information and that uses coherent-evolution-on-columns biclusters. This alternative approach shares some common features with sequential pattern mining. Finally, given a dataset of visitor trajectories, we indicate how these approaches can be used to build a collaborative recommendation strategy.

Introduction

CrossCult (http://www.crosscult.eu) is a European project whose idea is to support the emergence of a European cultural heritage by allowing visitors in different cultural sites (e.g. museum, historic city, archaeological site) to improve the quality of their visit by using adapted computer-based devices and to consider the visit at a European level. Such improvement can be accomplished by studying, among others, the possibility to build a dynamic recommendation system. This system should be able to produce a relevant suggestion on which part of a cultural site may be interesting for a specific visitor.

Here, our objective is to study a dynamic recommendation system for visitors in a museum. Given a new visitor V n , the task is to suggest a museum item that may be interesting for him/her. Based on how a suggestion is made to a new visitor V n , a recommendation system can be classified into one of the three following categories [START_REF] Adomavicius | Toward the next generation of recommender systems: A survey of the state-of-the-art and possible extensions[END_REF]:

-Content-based recommendations: The system makes a suggestion based only on the previous visited items of V n . For example, if V n visited mostly the items from prehistoric era, then the system recommends another item from that era.

-Collaborative recommendations: The system looks for previous users who have similar interest to V n , and makes a suggestion based on their visited items. For example, if many of V n 's similar users have visited item I, then the system recommends this item. -Hybrid approaches: The combination of content-based and collaborative approaches.

Our method belongs to the second category (collaborative recommendation).

First we group all previous users based on their visit trajectories using biclustering. When V n arrives, we try to find a G s , i.e. a group of visitors who shares a similar interest to V n . Then, based on the behavior of the visitors in G s , we can suggest one item that may be interesting for V n .

In this paper we will recall an approach in [7] that uses partition pattern structures to obtain biclusters with constant (or similar) values on the columns. Then we will propose an alternative approach that relies on this approach to mine another type of biclusters: those with coherent evolution on the columns (CEC biclusters). This bicluster type is useful when we are dealing with a dataset of trajectories where each trajectory corresponds to an ordered list of items. Furthermore, the mining of CEC biclusters can be related to sequential pattern mining, which we will explore in this paper.

This paper is organized as follows. First, we mention some related works about recommendation in Section 2. Then the basic background on biclustering is given in Section 3. Section 4 explains how to perform biclustering using partition pattern structures. The application of biclustering to recommendation systems will be presented in Section 5. Finally, we conclude our paper and outline some future works in Section 6.

Related work

In this section, we will describe related work about recommendation systems, biclustering, and Formal Concept Analysis (FCA).

FCA has been studied in collaborative movie recommendations for a user by looking at the ratings given by other users. In [START_REF] Du Boucher-Ryan | Collaborative recommending using formal concept analysis[END_REF], FCA is used to generate a lattice from a binary matrix (with users as rows and movies as columns) as the formal context. This matrix is derived from a rating dataset which is binarized, such that the matrix contains only the information whether a user has rated a movie. The lattice is then drawn to select some neighbors -i.e. users who have rated the same movies as the new user -regardless of the rating values. In this way, the exhaustive search of neighbors can be avoided. The neighbors' ratings can be then studied to recommend movies rated by the neighbors but not yet rated by the new user.

Pattern structures [START_REF] Ganter | Pattern structures and their projections[END_REF][START_REF] Kaytoue | Biclustering meets triadic concept analysis[END_REF] are a generalization of FCA, where the objects have more complex descriptions (e.g. sequence, graph, etc.). FCA was also extended into Triadic Concept Analysis, and it was shown in [START_REF] Kaytoue | Biclustering meets triadic concept analysis[END_REF] that triadic concepts are in 1-1-correspondence with maximal biclusters of similar values.

Partition pattern structures are an instance of the pattern structure framework. They were used in a collaborative movie recommendation [START_REF] Codocedo-Henríquez | Contributions à l'indexation et à la récupération d'information utilisant l'analyse formelle de concepts[END_REF] by identifying similar-column biclusters within the rating matrix. Such a bicluster corresponds to a set of users with similar rating behavior (hence similar interest) across a set of movies. Therefore, to recommend a movie to a new user, there is a search for biclusters whose users have similar interest to him/her. Using the real MovieLens data, a study was also conducted based on Boolean matrix factorization [START_REF] Akhmatnurov | Context-aware recommender system based on boolean matrix factorisation[END_REF].

Moreover, recommendation systems based on FCA and/or biclustering have been applied to other real world problems such as detection of future advertising terms for a company [START_REF] Ignatov | Concept-Based Biclustering for Internet Advertisement[END_REF], educational orientation of Russian school graduates [START_REF] Ignatov | Recommender system based on algorithm of bicluster analysis RecBi[END_REF], and idea recommendation at a crowdsourcing project of Witology company [START_REF] Ignatov | FCA-based recommender models and data analysis for crowdsourcing platform Witology[END_REF]. Other than recommendation systems, biclustering has also been applied in other fields, for example in the study of miRNA-gene target interaction data and miRNA functions and mechanisms [START_REF] Pio | ComiRNet: a web-based system for the analysis of miRNA-gene regulatory networks[END_REF]. A unified taxonomy of biclustering methods was proposed in [START_REF] Ignatov | Towards a unified taxonomy of biclustering methods[END_REF].

Biclustering

In this section, we will recall the basic background and discuss illustrative examples of the different types of biclusters as described in [START_REF] Madeira | Biclustering algorithms for biological data analysis: a survey[END_REF].

We consider a dataset composed of a set of objects, each of which has values over a set of attributes. This dataset can be represented as a numerical matrix, where each cell ij indicates the value of object i w.r.t. attribute j.

One may be interested in finding which subset of objects possesses the same values w.r.t. a subset of attributes. Regarding the matrix representation, this is equivalent to the problem of finding a submatrix that has a constant value over all of its elements (example in Table 1). This task is called biclustering with constant values, which is a simultaneous clustering of the rows and columns of a matrix. Other than constant values, the bicluster approach also focused on finding other types of submatrices, as shown in Table 2. A bicluster with constant columns (rows) is a submatrix where each column (row) has the same value, as illustrated in Table 2a (Table 2b, resp.).

In a bicluster with additive coherent values, the value of each cell ij follows the equation γ + α i + β j , where γ is a constant, α i is a constant value for row i, and β j is a constant value for column j. For example, if γ = 1, (α 1 , α 2 , α 3 , α 4 ) = [START_REF] Alqadah | Biclustering neighborhood-based collaborative filtering method for top-n recommender systems[END_REF][START_REF] Akhmatnurov | Context-aware recommender system based on boolean matrix factorisation[END_REF][START_REF] Ben-Dor | Discovering local structure in gene expression data: the order-preserving submatrix problem[END_REF][START_REF] Codocedo | Lattice-based biclustering using partition pattern structures[END_REF], and (β 1 , β 2 , β 3 , β 4 ) = (0, -2, 1, -1), then we can obtain the bicluster in Table 2c. Similarly, we can obtain a bicluster with multiplicative coherent values as shown in Table 2d using a constant for each row and each column. The main difference is that, instead of adding, we multiply them. Another interesting type is the CEC bicluster, also known as order-preserving submatrix [START_REF] Ben-Dor | Discovering local structure in gene expression data: the order-preserving submatrix problem[END_REF]. In this type of bicluster, each row induces the same linear order across all columns. For example, in the bicluster in Table 2e, each row follows column1 ≤ column2 ≤ column4 ≤ column3. Moreover, a bicluster with coherent evolution on the rows can be defined similarly, as shown in Table 2f.

Those different types of biclusters are useful when we are interested in identifying a group of people who behave similarly according to a set of attributes. This group identification is necessary in the task of collaborative recommendation, because in the process of making a suggestion to a person, we first identify the people who are similar to him/her.

FCA and Pattern Structures

(FCA) is a mathematical framework based on lattice theory and used for classification, data analysis, and knowledge discovery [START_REF] Ganter | Formal Concept Analysis: Mathematical Foundations[END_REF]. From a formal context, FCA detects all formal concepts, and arranges them in a concept lattice. Definition 1. A formal context is a triple (G, M, I), where G is a set of objects, M is a set of attributes, and I is a binary relation between G and M , i.e.

I ⊆ G × M .
If an object g has an attribute m, then (g, m) ∈ I. An example of a formal context is shown in Table 3. This table shows whether a visitor (V 1 -V 4 ) visits an item (102, 302, 402, or 704).

The Galois connection for a formal context (G, M, I) is defined as follows:

Definition 2. For a subset of objects A ⊆ G, A is the set of attributes that are possessed by all objects in A, i.e.:

A = {m ∈ M |∀g ∈ A, (g, m) ∈ I}, A ⊆ G Dually,
V1 × V2 × × V3 × × × V4 × × × A formal concept (A, B) is a subconcept of (C, D) -denoted by (A, B) ≤ (C, D) -if A ⊆ C (or equivalently D ⊆ B).
A concept lattice can be formed using the ≤ relation which defines the order among concepts. For the context in Table 3, the formal concepts and their corresponding lattice are shown in Fig. 1. 3 FCA is restricted to specific datasets where each attribute is binary (e.g. has only yes/no value). For more complex values (e.g. numbers, strings, trees, graphs. . . ), FCA is then generalized into pattern structures [START_REF] Ganter | Pattern structures and their projections[END_REF]. The operator is a similarity operation that returns the common elements between any two descriptions. A description can be a set, a sequence, or other complex structure. In the case of set as a description, corresponds to set intersection (∩), i.e. {a, b, c} {a, b, d} = {a, b}, and corresponds to subset inclusion (⊆). In the case of sequence as a description, is a set of common closed subsequences (SCCS) [START_REF] Codocedo | Lattice-based biclustering using partition pattern structures[END_REF]. Similarly, corresponds to subsequence inclusion ( ).

Biclustering Using Partition Pattern Structures

Biclustering has many common elements with FCA. In FCA, from a binary matrix we try to find a maximal submatrix whose elements are 1. In other words, the objective is to identify maximal constant-value biclusters (but only for biclusters whose values are 1). Hence, a formal concept can be considered as a bicluster of objects and attributes. Furthermore, formal concepts are arranged in a concept lattice, that can describe the hierarchical relation among all biclusters. Consider the matrix given by Table 4, where we are interested in finding constant-column biclusters. We recognize that the values of m 1 "break" the objects into two sets: {g 1 , g 2 } and {g 3 , g 4 }. The same "break" is also obtained from the values of m 4 . In particular, we can see that the pair ({g 1 , g 2 }, {m 1 , m 4 }) corresponds to a constant-column bicluster. Therefore, it is possible to mine this type of bicluster using this "breaking" -or "partitioning" -technique. Moreover, this technique can be performed using partition pattern structures -an extension of FCA.

In this section, first we will recall the constant-column biclustering approach using partition pattern structures [7]. We will then propose an extension of this approach to perform CEC biclustering.

Biclustering with Constant Columns

A partition d = {p i } of a set G is a collection of p i ⊆ G such that: pi∈d p i = G and p i ∩ p j = ∅ whenever i = j. (1)
Notice that when calculating the initial partitions, missing values can produce overlapping partitions (i.e. p i ∩ p j = ∅). Consider the dataset given by Table 4 that has G = {g 1 , g 2 , g 3 , g 4 } as the set of objects and M = {m 1 , m 2 , m 3 , m 4 , m 5 } as the set of attributes. Here we can define a partition mapping δ : M → D. The partition is based on the fact that the values of the attribute are equal for all objects in a subset. For example, δ(m 1 ) = {{g 1 , g 2 }, {g 3 , g 4 }} because G is partitioned as such regarding the value of m 1 , whereas δ(m 4 ) = {{g 1 , g 2 }, {g 1 , g 3 , g 4 }}. This partition overlaps on g 1 since this object has a missing value on m 4 . Intuitively, g 1 can be grouped with either {g 2 } or {g 3 , g 4 } w.r.t. m 4 .

The space D of all partitions over G is a complete lattice, where the meet and join of two partitions d 1 = {p i } and d 2 = {p j } are defined as: Table 4. A dataset with 4 objects and 5 attributes m1 m2 m3 m4 m5 g 1 1 5 3 ? 7 g 2 1 1 4 2 7 g 3 2 5 4 5 3 g 4 2 5 4 5 7

d 1 d 2 =   i,j p i ∩ p j   + (2)
d 1 d 2 =   pi∩pj =∅ p i ∪ p j   + (3)
where (.) + is a closure that preserves only the maximal components in d. For example, δ(m

1 ) δ(m 4 ) = {{g 1 , g 2 }, {g 1 }, {g 3 , g 4 }} + = {{g 1 , g 2 }, {g 3 , g 4 }}, and δ(m 1 ) δ(m 4 ) = {{g 1 , g 2 }, {g 1 , g 2 , g 3 , g 4 }, {g 1 , g 3 , g 4 }} + = {{g 1 , g 2 , g 3 , g 4 }}.
The order between any two partitions is given by the subsumption relation:

d 1 d 2 ⇐⇒ d 1 d 2 = d 1 (4)
Given a set of attributes M, a partition space D, and a mapping δ, a partition pattern structures for constant-column biclustering is determined by the triple (M, (D, ), δ). A pair (A, d) is then called a partition pattern concept (pp-concept) iff A = d and d = A, where:

A = m∈A δ(m) A ⊆ M (5)
d = {m ∈ M|d δ(m)} d ∈ D (6)
For any partition component p ∈ d, each pair (p, A) corresponds to a constantcolumn bicluster. For example, from the concept ({m 1 , m 4 }, {{g 1 , g 2 }, {g 3 , g 4 }}), two biclusters can be obtained: ({g 1 , g 2 }, {m 1 , m 4 }) and ({g 3 , g 4 }, {m 1 , m 4 }).

Biclustering with Coherent Evolution on the Columns

In a dataset of movie ratings, bicluster with constant columns is useful to identify a set of users with the same taste regarding a set of movies. Another interesting problem arises, e.g. when the dataset contains watching order. In that case, we may be interested in finding a set of users who watch a set of movies in the same order. This problem corresponds to CEC biclustering, where the objective is to find a set of rows which has coherent evolution over a set of columns, as Table 5. A dataset with 5 objects and 5 attributes m1 m2 m3 m4 m5 g 1 1 2 3 4 5 g 2 4 2 1 ? 3 g 3 2 3 4 1 1 g 4 5 4 2 3 1 g 5 2 1 5 4 3 Table 6. Some examples of partitions over Table 5 Pair Partition p 1,2 {{g 1 , g 3 }, {g 2 , g 4 , g 5 }} p 1,3 {{g 1 , g 3 , g 5 }, {g 2 , g 4 }} p 1,4 {{g 1 , g 2 , g 5 }, {g 2 , g 3 , g 4 }} p 2,3 {{g 1 , g 3 , g 5 }, {g 2 , g 4 }} p 2,5 {{g 1 , g 2 , g 5 }, {g 3 , g 4 }} previously described in Section 3. In the current section, we will explain the possible application of partition pattern structures to discover CEC biclusters.

Consider the dataset given by Table 5, with the set of attributes G = {g 1 , g 2 , g 3 , g 4 , g 5 }. First, we have to list each pair of attributes and the partition according to the pair's evolution. For the pair p 1,2 = (m 1 , m 2 ), the partition is {{g 1 , g 3 }, {g 2 , g 4 , g 5 }} because in g 1 and g 3 , m 1 is less than m 2 , whereas in g 2 , g 4 , and g 5 , m 1 is greater. As in Subsection 4.1, missing values generate an overlapping partition. For instance, the pair p 1,4 gives rise to the partition is {{g 1 , g 2 , g 5 }, {g 2 , g 3 , g 4 }}. Furthermore, two columns with the same value (e.g. g 3 in m 4 and m 5 ) can also produce an overlapping partition because, by our definition of CEC bicluster in Section 3, they satisfy m 4 ≤ m 5 and m 5 ≤ m 4 . Therefore, the partition for p 4,5 is {{g 1 , g 2 , g 3 }, {g 2 , g 3 , g 4 , g 5 }}. Some pairs and their partitions are listed in Table 6.

Since a CEC partition is defined by at least two attributes, the partition mapping becomes γ : P → D. For instance, γ(p 1,2 ) = {{g 1 , g 3 }, {g 2 , g 4 , g 5 }}.

As in Subsection 4.1, given a set of attribute pairs P, a partition space D, and the mapping function γ, a partition pattern structures for coherent-evolution biclustering is determined by the triple (P, (D, ), γ). A pp-concept is a pair (B, d) such that B = d and d = B, where:

B = p∈B γ(p) B ⊆ P (7) d = {p ∈ P|d γ(p)} d ∈ D (8)
Here, the extent of a pp-concept is a set of attribute pairs. We can obtain a CEC bicluster in a pp-concept if there is a clique among the attributes in the pairs. For example, consider the pp-concept ppc 1 with extent {p 1,2 , p 1,3 , p 2,3 } and intent {{g 1 , g 3 }, {g 5 }, {g 2 , g 4 }}. Its extent forms a clique among m 1 , m 2 , and m 3 , since all pairings of any two of those attributes are included.

If a pp-concept (B, d) contains a set of attributes A that forms a clique, then each pair (p, A), for any partition component p ∈ d, corresponds to a CEC bicluster. For example, from ppc 1 , we can obtain bicluster ({g 1 , g 3 }, {m 1 , m 2 , m 3 }).

Comparison with Sequential Pattern Mining

A sequence is an ordered list s 1 s 2 . . . s m , where s i is an itemset {i 1 , . . . , i n }. A sequence s = s 1 s 2 . . . s m is a subsequence of s = s 1 s 2 . . . s n , denoted by s s , if there exist indices 1 ≤ i 1 < i 2 < . . . < i m ≤ n such that s j ⊆ s ij for all j = 1 . . . m and m ≤ n. For example, the sequence {a}{d} is a subsequence of {a, b}{a, c, d} , while sequence {c}{d} is not.

Notice that the problem of retrieving CEC biclusters can be thought of as a particular type of sequential pattern mining where each itemset is composed by only one item, i.e. the sequences are an ordered list of items. Mining sequential patterns means retrieving frequent subsequences (i.e., subsequences that are present in more than n sequences) and for which there exist many efficient algorithms [START_REF] Yan | CloSpan: Mining: Closed sequential patterns in large datasets[END_REF][START_REF] Wang | BIDE: Efficient mining of frequent closed sequences[END_REF].

The CEC biclustering differs from sequential pattern mining when we allow overlaps in the partitions. Consider Table 5 as a sequential dataset. Each number in row x column y corresponds to the itemset when item y appears in the sequence x. For example, the sequence of object g 3 is {m 4 , m 5 }{m 1 }{m 2 }{m 3 } .

Let us consider items m 1 , m 4 , and m 5 . According to sequential pattern mining, g 3 is different from g 4 , because m 4 and m 5 appear in the same itemset in g 3 . On the other hand, according to CEC biclustering with overlaps, g 3 is similar to g 4 , because in both objects m 5 ≤ m 4 ≤ m 1 .

Recommendation

In the context of CrossCult, we are working on a visitor dataset that comprises several trajectories in a museum. Within this project, our main objective is to build a dynamic recommendation system for new visitors. This system should be able to suggest a museum item to visitors based on their trajectories and by looking at the trajectories of previous visitors. Also, it should be able to update the suggestion as they move inside the museum.

Matrix as order of interest

For each item in the museum, we can measure (e.g. by rating, duration of visit, etc.) their level of interestingness from a set of visitors. An example is shown in Table 7, where the number in cell xy is the ranking of item y according to visitor x. Here we have 3 visitors (v 1 , v 2 , and v 3 ) in the database and 1 target visitor (v a ). Among the existing visitors, only v 1 has complete values over all five items. According to this visitor i 1 is the best, followed by i 2 , i 3 , i 4 , and the worst i 5 . For v 2 (v 3 ), the order of interest of i 2 (i 4 resp.) is not known.

The target visitor (v a ) has visited only three items, with the same order of preference in i 1 and i 2 . This visitor will be included in the bicluster ({v 1 , v 2 , v a }, {i 1 , i 2 , i 3 }). The other two members of this bicluster do not agree on the order of interest of i 4 and i 5 . Hence, we should suggest i 5 to him/her, since one visitor (v 2 ) similar to him/her ranked it first.

Matrix as order of visit

Consider the dataset given by Table 8 dataset. The numbers in row x indicate the path of visitor v x . For example, the path of

v 2 is i 6 → i 1 → i 4 → i 2 → i 3 → i 7 → i 5 .
Now we have a new visitor v a who recently arrives to the museum. He/She visits i 3 , followed by i 4 . Our task is to recommend a new item to her, by studying the CEC biclusters over the first four visitors. Some of those biclusters are listed in Table 9. From B7, we can see that all four visitors visit i 7 after i 2 . In B3, v 1 and v 3 follow the same order w.r.t. {i 3 , i 4 , i 5 , i 7 }: i 3 → i 4 → i 5 → i 7 . These two visitors also agree in the order of items {i 1 , i 4 , i 5 , i 7 }, as seen in B4. Those CEC biclusters can be studied to give a recommendation to v a by focusing on those visitors that are similar to him/her. Thus, we can propose a recommendation strategy that follows sequential patterns in the dataset. The idea behind is the following: if many visitors have the path i a → i b → i c , then we should recommend item i c to a visitor who has done i a → i b .

Since v a has path i 3 → i 4 , we focus on the CEC biclusters that have those two items, i.e. B1, B2, and B3. One of those biclusters (B2) has a different ordering (i 3 after i 4 ), and thus we filter it out. Then, in B3 for example, the path is i 3 → i 4 → i 5 → i 7 . Therefore, we can recommend i 5 to v a .

Application to Hecht Museum

In the framework of the CrossCult project, we are working on a specific dataset about the trajectories of 254 visitors in Hecht Museum in Haifa, Israel [START_REF] Lanir | The influence of a locationaware mobile guide on museum visitors' behavior[END_REF]. In this dataset, there are 52 items over 8 rooms (A-G). Therefore, the visitor-item matrix for biclustering will have 254 rows and 52 columns.

A visitor's trajectory corresponds to a list of visits, which is composed by three elements: "start time", "end time", and "item name". From the trajectories, we can build either matrix of order of interest or matrix of order of visit. Then, the application of recommendation can be tested over this dataset.

Conclusion

In this work, we have explored an approach to build collaborative recommendation strategy for visitors in a museum. This strategy takes into account the order of interest or the order of visit for each visitor, and we showed how to use CEC biclustering to obtain a set of similar visitors. We also presented a technique for mining CEC biclusters based on FCA using partition pattern structures. This recommendation strategy can be applied to any dataset where the order of items is relevant.

As future work, we intend to explore recommendations based on the order of visit. The problem to be solved is how to model visitors who visit a single item multiple times (for example, i 1 , i 2 , and back to i 1 ).

Another noteworthy question is how to measure the "score" of each bicluster in order to rank recommendations for a new visitor. Ranking candidate items was studied in [START_REF] Alqadah | Biclustering neighborhood-based collaborative filtering method for top-n recommender systems[END_REF] for constant-value biclusters, and it is possible to extend this work to CEC biclusters. Moreover, further comparisons of CEC biclustering and sequential pattern mining should be investigated, in particular, regarding their complexities and their results. Biclustering-based recommender systems could also be compared to social network analysis and cluster-indexing collaborative filtering [START_REF] Kim | Recommender systems using cluster-indexing collaborative filtering and social data analytics[END_REF]. Finally, an implementation of the CEC biclustering using partition pattern structures and an empirical study on real-world data should be performed to measure its complexity and efficiency.

Table 2 .

 2 Examples of some types of biclusters. (a) Constant columns, (b) constant rows, (c) additive coherent values, (d) multiplicative coherent values, (e) coherent evolution on the columns, and (f) coherent evolution on the rows.

  for a subset of attributes B ⊆ M , B is the set of objects that have all attributes in B, i.e.: B = {g ∈ G|∀m ∈ B, (g, m) ∈ I}, B ⊆ M Definition 3. A formal concept is a pair (A, B), where A ⊆ G and B ⊆ M , and such that A = B and B = A.

Fig. 1 .

 1 Fig. 1. Concept lattice for the formal context in Table3

Definition 4 .

 4 A pattern structure is a triple (G, (D, ), δ), where G is a set of objects, (D, ) is a complete meet-semilattice of descriptions, and δ : G → D maps an object to a description.

Definition 5 .

 5 The Galois connection for a pattern structure (G, (D, ), δ) is defined as:A = g∈A δ(g), A ⊆ G d = {g ∈ G|d δ(g)}, d ∈ DFinally, a pattern concept is similar to a standard formal concept: Definition 6. A pattern concept is a pair (A, d), A ⊆ G and d ∈ D, where A = d and d = A.

Table 1 .

 1 A

	bicluster with constant value (shaded)
	1 1 4 3 5
	1 1 2 5 1
	3 3 4 2 1

Table 3 .

 3 A formal context for four objects, with four items: 102, 302, 402, and 704 as an example.

	102 302 402 704

  about 4 visitors in a museum with 7 items. As explained in Subsection 4.3, this table can be regarded as a sequential

Table 7 .

 7 Order of interest of 5 items, observed from certain visitors

	i1 i2 i3 i4 i5
	v1 1 2 3 4 5
	v2 3 ? 1 2 4
	v3 2 4 3 ? 1
	va 1 1 2 ? ?

Table 8 .

 8 Order of visit of 7 items, observed from certain visitors

	i1 i2 i3 i4 i5 i6 i7
	v1 1 2 3 4 5 6 7
	v2 2 4 5 3 7 1 6
	v3 4 2 1 5 6 3 7
	v4 7 3 1 4 2 6 5
	va ? ? 1 2 ? ? ?

Table 9 .

 9 Some CEC biclusters in Table8

	# Visitors	Items (in order)
	B1 v1	i1, i2, i3, i4, i5, i6, i7
	B2 v2	i6, i1, i4, i2, i3, i7, i5
	B3 v1, v3	i3, i4, i5, i7
	B4 v1, v3	i1, i4, i5, i7
	B5 v1, v4	i3, i5, i6
	B6 v2, v3	i6, i1, i4, i5
	B7 v1, v2, v3, v4 i2, i7