[PDF] A smoothed and probabilistic PARAFAC model with covariates

View - Download A smoothed and probabilistic PARAFAC model with covariates

Previous PDF

Next PDF

Smooth and probabilistic PARAFAC model

with auxiliary covariates

Leying Guan

Department of Biostatistics, Yale University

August 30, 2022

Abstract

In immunological and clinical studies, matrix-valued time-series data clustering is increasingly popular. Researchers are interested in nding low-dimensional embed- ding of subjects based on potentially high-dimensional longitudinal features and in- vestigating relationships between static clinical covariates and the embedding. These studies are often challenging due to high dimensionality, as well as the sparse and ir- regular nature of sample collection along the time dimension. We propose a smoothed probabilistic PARAFAC model with covariates (SPACO) to tackle these two prob- lems while utilizing auxiliary covariates of interest. We provide intensive simulations to test dierent aspects of SPACO and demonstrate its use on an immunological data set from patients with SARs-CoV-2 infection. Keywords:Tensor decomposition; Time series; Missing data; Probabilistic model.

1arXiv:2?04.05?84v3 [stat.ME] 27 Aug 2022

1In troduction

Sparsely observed multivariate times serious data are now common in immunological stud- ies. For each subject or participant i (i= 1;:::;I), we can collect multiple measurements onJfeatures over time, but often atnidierent time stampsfti;1;:::;ti;nig. For exam- ple, for each subject, immune proles are measured for hundreds of markers at irregular sampling times in Lucas et al. (2020) and Rendeiro et al. (2020). LetXi2RniJbe the longitudinal measurements for subjecti, we can collectXifor allIsubjects into a sparse three-way tensorX2RITJ, whereT=j [ifti;1;:::;ti;nigjis the number of distinct time stamps across all subjects. Sincefti;1;:::;ti;nigtend to be small in size and have low overlaps for dierent subjecti,Xmay have a high missing rate along the time dimension. In addition toXi, researchers often have a set of nontemporal covariateszi2Rqfor subjectisuch as medical conditions and demographics, which may account partially for the variation in the temporal measurementsXacross subjects. Modeling such auxiliary covariatesZ:= (z1;:::;zI)>together withXmight help with the estimation quality and understanding for the cross-subject heterogeneity. In this paper, we propose SPACO (smooth and probabilistic PARAFAC model with auxiliary covariates) to adapt to the sparsity long the time dimension inXand utilize the auxiliary variablesZ. SPACO assumes thatXis a noisy realization of some low-rank signal tensor whose time components are smooth and subject scores have a potential dependence on the auxiliary covariatesZ: x itj=KX k=1u iktkvjk+itj; itj N(0;2j) u i= (uik)Kk=1 N(i;f);i=>zi: Here, (1)2RqKdescribes the dependence of the expected subject scoreifor subject ionzi, and (2)uik,tk,vjkare the subject score, trajectory value and feature loading for factorkin the PARAFAC model and the observation indexed by (i;t;j) whereuihas a normal priorN(i;f). We impose smoothness on time trajectories (tk)>t=1and sparsity onto deal with the irregular sampling along the time dimension and potentially high dimensionality inZ. 2 Alongside the model proposal, we will also discuss several issues related to SPACO, including model initialization, auto-tuning of smoothness and sparsity inand hypothesis testing onthrough cross-t. Successfully addressing these issues is crucial to practitioners interested in applying SPACO to their analysis. In the remaining of the article, we summarize some closely related work in section 1.1 and describe the SPACO model in Section 2 and model parameter estimation with xed tuning parameters in Section 3. In Section 4, we discuss the aforementioned related issues that could be important in practice. We compare SPACO to several existing methods in Section 5. Finally, in Section 6, we apply SPACO to a highly sparse tensor data set on immunological measurements for SARS-COV-2 infected patients. We provide a python packageSPACOfor researchers interested in applying the proposed method. 1.1

Related w ork

In the study of multivariate longitudinal data in economics, researchers have combined ten- sor decomposition with auto-cross-covariance estimation and autoregressive models (Fan et al., 2008; Lam et al., 2011; Fan et al., 2011; Bai and Wang, 2016; Wang et al., 2019,

2021). These approaches do not work well with highly sparse data or do not scale well

with the feature dimensions, which are important for working with medical data. Func- tional PCA (Besse and Ramsay, 1986; Yao et al., 2005) is often used for modeling sparse longitudinal data in the matrix-form. It utilizes the smoothness of time trajectories to handle sparsity in the longitudinal observations, and estimates the eigenvectors and factor scores under this smoothness assumption. Yokota et al. (2016) and Imaizumi and Hayashi (2017) introduced smoothness to tensor decomposition models, and estimated the model parameters by iteratively solving penalized regression problems. The methods above don't consider the auxiliary covariatesZ. It has been previously discussed that includingZcould potentially improve our esti- mation. Li et al. (2016) proposed SupSFPC (supervised sparse and functional principal component) and observed that the auxiliary covariates improve the signal estimation qual- ity in the matrix setting for modeling multivariate longitudinal observations. Lock and Li (2018) proposed SupCP which performs supervised multiway factorization model with 3 complete observation and employs a probabilistic tensor model (Tipping and Bishop, 1999; Mnih and Salakhutdinov, 2007; Hinrich and Mrup, 2019). Although an extension to sparse tensor is straightforward, SupCP does not model the smoothness and can be much more aected by severe missingness along the time dimension. SPACO can be viewed as an extension of functional PCA and SupSFPC to the setting of three-way tensor decomposition (Acar and Yener, 2008; Sidiropoulos et al., 2017) using the parallel factor (PARAFAC) model (Harshman and Lundy, 1994; Carroll et al., 1980). It uses a probabilistic model and jointly models the smooth longitudinal data with potentially high-dimensional non-temporal covariatesZ. We refer to the SPACO model as SPACO- when no auxiliary covariateZis available. SPACO- itself is an attractive alternative to existing tensor decomposition implementations with probabilistic modeling, smoothness regularization, and automatic parameter tuning. In our simulations, we compare SPACO with SPACO to demonstrate the gain from utilizingZ. 2

SP ACOMo del

2.1

Notations

LetX2RITJbe a tensor for some sparse multivariate longitudinal observations, where Iis the number of subjects,Jis the number of features, andTis the number of total unique time points. For any matrixA, we letAi:=A:idenote itsithrow/column, and often write A :iasAifor theithcolumn for convenience. LetXI:= X :;:;1X:;:;J

2RI(TJ),

X T:= X >:;:;1X>:;:;J

2RT(IJ),XJ:=

X >:;1;:X>:;T;:

2RJ(IT)be the matrix

unfolding ofXin the subject/feature/time dimension respectively. We also dene: Tensor product}:a2RI,b2RT,c2RJ, then,A=a}b}c2RITJwith A itj=aibtcj.

Kronecker product

:A2RI1K1,B2RI2K2, then C=A B=0 B BB@A

11B ::: A1K1B

A

I11B ::: AI1K1B1

C

CCA2R(I1I2)(K1K2):

Column-wise Khatri-Rao product:A2RI1K,B2RI2K, thenC=AB2 4 R (I1I2)KwithC:;k= (A:;k

B:;k) fork= 1;:::;K.

Element-wise multiplication:A;B2RIK, thenC=AB2RIKwithCik= (AikBik); forb2RK,C=Ab=Adiagfb1;:::;bKg; forb2RI,C=bA= diagfb1;:::;bIgA. 2.2 smo othand probabilistic P ARAFACmo delwith co variates We assumeXto be a noisy realization of an underlying signal arrayF=PK k=1Uk}k}Vk. We stackUk=k=Vkas the columns ofU==V, denote the rows ofU==Vbyui=t=vj, and their entries byuik=tk=vjk. We letxitjdenote the (i;t;j)-entry ofX. Then, x itj=KX k=1u iktkvjk+itj;ui N(i;f); ijt N(0;2j);(1) where f= diagfs21;:::;s2Kgis aKKdiagonal covariance matrix. Even thoughXis often of high rank, we consider the scenario where the rankKofFis small. Without covariates, we set the prior mean parameteri=0. If we are interested in explaining the heterogeneity iniacross subjects with auxiliary covariatesZ2RIq, then we may modelias a function ofzi:=Zi;:. Here, we consider a linear model ik=z>ik;8k= 1;:::;K. To avoid confusion, we will always callXthe \features", andZthe \covariates" or \variables". We refer toUas the subject scores, which characterize dierences across subjects and are latent variables. We refer toVas the feature loadings, which reveal the composition of the factors using the original features and could assist the downstream interpretation. Finally,is referred to as the time trajectories, which can be interpreted as function values sampled from some underlying smooth functionsk(t) at a set of discrete-time points, e.g., k= (k(t1);:::;k(tT)). Recalling thatXI2RI(TJ)is the unfolding ofXin the subject direction, we write ifor the indices of observed values in theithrow ofXI, andXI;~ifor the vector of these observed values. Each such observed valuexitjhas noise variance2j, and we write ~i to represent the diagonal covariance matrix with diagonal values2jbeing the corre- sponding variances for"ijtat indices in~i. Similarly, we denef~t;XT;~t;~tgfor the unfoldingXT2RT(IJ), andf~j;XJ;~j;~jgfor the observed indices, the associated ob- 5 served vector and diagonal covariance matrix for thejthrow inXJ2RJ(IT). We set =fV;;;2j;j= 1;:::;J;(s2k;k= 1;:::;K)gto denote all model parameters. Set H= (V) andfi=XI;~iH~ii. IfUis observed, the complete data log-likelihood is

L(X;Uj) =12

X i f >i1 ~ifi+ (uii)>1 f(uii) + logj~ij+Ilogjfj :(2) Set ~i= ~i+H~ifH>~i. The marginalized log likelihood integrating out the randomness inUenjoys a closed form (Lock and Li, 2018):

L(Xj)/ 12

X if >i~1 ifi+ logj~ij! :(3) Set i= (H>~i1 ~iH~i+ 1 f)1. We can also equivalent express the marginal likelihood as below.

L(Xj)/ 12

f>i 1 ~i1 ~iH~iiH>~i1 ~i f i12 (logj~ij+ logjfj logjij);(4) We use the form in eq. (4) to derive the updating formulas and criteria for rank selection since it does not involve the inverse of a large non-diagonal matrix. Model parameters in eq. (3) or eq. (4) are not identiable due to (1) parameters rescaling from (k;Vk;k;s2k) to (c1k;c2Vk;c3k;c23s2k) for anyc1c2c3= 1, and (2) reordering of dierent componentkfork= 1;:::;K. More discussions of the model identiability can be found in Lock and Li (2018). Hence, adopting similar rules from Lock and Li (2018), we require (C.1)kVkk22= 1,kkk22=T. (C.2) The laten tcomp onentsare in decreasing order based on their o verallv ariances f;kk+ kZkk22=I, and the rst non-zero entries inVkandkto be positive, e.g.,vk1>0 andk1>0 if they are non-zero. To deal with the sparse sampling along the time dimension and take into consideration that features are often smooth over time in practice, we assume that the time component kis sampled from a slowly varying trajectory functionk(t), and encourage smoothness ofk(t) by directly penalizing the function values via a penalty termP k1k>k k. This 6 paper considers a Laplacian smoothing (Sorkine et al., 2004) with a weighted adjacency matrix . LetT(t) represent the associated time forX:;t:;. We dene and as >; =0 B

BBBBB@1T[2]T[1]1T[2]T[1]:::0 0

0

1T[3]T[2]...0...

0 0:::1T[T]T[T1]1T[T]T[T1]1

C

CCCCCA2RT(T1)

Practitioners may choose other forms for

. If practitioners wantk(t) to have slowly varying derivatives, they can also use a penalty matrix that penalizes changes in gradients over time. Further, when the number of covariatesqinZis moderately large, we may wish to impose sparsity in theparameter. We encourage such sparsity by including a lasso penalty (Tibshirani, 2011) in the model. In summary, our goal is then to nd parameters maximizing the expected penalized log-likelihood, or minimizing the penalized expected deviance loss, under norm constraints: minJ():=12

L(Xj) +KX

k=1 1k>k k+X k 2kjkj s.tkVkk22= 1;kkk22=T;for allk= 1;:::;K:(5) Only the identiability constraint (C.1) has entered the objective. We can always guarantee (C.2) by changing the signs inV,,and reordering the components afterward without changing the achieved objective value. Eq. (5) describes a non-convex problem. We will nd locally optimal solutions via an alternating update procedure: (1) xing other parameters and updatingvia lasso regressions; (2) xingand updating other model parameters using the EM algorithm. We give details of our iterative estimation procedure in Section 3. 3

Mo delparameter estimation

Given the model rankKand penalty terms1k,2k, we alternately update parameters, V,U,s2and2with a mixed EM procedure described in Algorithm 1. We brie y explain the updating steps here: 7 (1):Given other parameters, we ndto to directly minimize the objective by solving a least-squares regression problem with lasso penalty. (2):Fixing, we update the other parameters using an EM procedure. Denote the current parameters as

0. Our goal is to minimize the penalized expected negative log-likelihood

J(;0):=EUj0(

L(X;Uj) +X

k 1k>k k+X k

2kjkj)

;(6) under the current posterior distributionUj0. We adopt a block-wise parameter updating scheme where we updateVk,k, fand2jsequentially.Algorithm 1:SPACO with xed penaltiesData:X, ,1,2,K Result:EstimatedV,,,s2,2and posteriorP(Uj;X) and the marginalized densityP(Xj) .

1Initialization ofV,,,s2,2and the posterior distribution ofU;

2whileNot convergeddo3fork= 1;:::;Kdo4:;k argmin:;kfL(Xj) +2kj:;kjg

5Vk argminVk[J(;0) +kVkk22],is the largest value leading to the

minimizer havingkVkk22= 1

6k argminVk[J(;0) +kUkk22],is the largest value leading to the

minimizer havingkkk22=T

7s2k argmins2kEUj0L(X;Uj).

8end

9Forj= 1;:::;J:2j argmin2jEUj0L(X;Uj).

10Update the posterior distribution ofU.

11endAlgorithm 1 describes the high level ideas of our updating schemes. In line 5 and

6, we guarantee the norm constraints onVkandkby adding an additional quadratic

term and set the coecientto guarantee the norm requirements. Even though this is not a convex problem, the proposed approaches provide optimal solutions for sub-routines updating dierent parameter blocks, and the penalized (marginalized) deviance loss is non- 8 increasing over the iterations. Theorem 3.1In Algorithm 1, let`and`+1are the estimated parameters at the begin- ning and end of the`thiteration of the outerwhileloop. We haveJ(`)J(`+1). Proof of Theorem 3.1 is given in Appendix B.1. In Algorithm 1, the posterior distribution ofuifor each row inUis Gaussian, with posterior covariance the same as idened earlier, and posterior mean given below. i= i 1 f>zi+ (H)>~i1 ~iXI;~i :(7) Explicit formulas and steps for carrying out the subroutines at lines 4-7 and line 9 are deferred to Appendix A.1 4

Initialization, tuning and testing

4.1

Initialization

One initialization approach is to form a Tucker decomposition [U?;?;V?;G] ofXusing HOSVD/MLSVD (De Lathauwer et al., 2000) whereG2RK1K2K3is the core tensor andU?2RIK1,?2RTK2,V?2RJK3are unitary matrices multiplied with the core tensors along the subject, time and feature directions respectively (K1=K2=K3is the smallest betweenKandI=T=J), and then perform PARAFAC decomposition on the small core tensorG(Bro and Andersson, 1998; Phan et al., 2013). We initialize SPACO with Algorithm 2, which combines the above approach with functional PCA (Yao et al., 2005) to work with sparse longitudinal data. Algorithm 2 consists of the following steps: (1) perform SVD onXJto getV?; (2) projectXJonto each column ofV?and perform functional PCA to estimate?; (3) run a ridge-penalized regression of rows ofXIonV? ?, and estimateU?andGfrom the regression coecients In a noiseless model with= 0 and complete temporal observations, one may replace the functional PCA step of Algorithm 2 with standard PCA. Then [U;;V] becomes a

PARAFAC decomposition of

11+X. 9

Algorithm 2:Initialization of SPACO1LetV?be the topK3left singular vectors ofXJusing only the observed columns.

2SetY(k) = (Y1(k);:::;YT(k))2RIT, whereYt(k) =X:;t;:(V?)k2RI.

3Let?be the topK2eigenvectors from functional PCA on the row aggregation of

matricesY(k)k= 1;:::;K3. (see Appendix A.2 for details on functional PCA.)

4Let~U= argminUfkXIU(V?

?)>k2F+kUk2Fg 2RIK2for some small.

5LetU?be the topKleft singular eigenvectors of~U, and~G=U>?~U2RKK2.

LetG2RKKKbe the estimated core array from rearranging~G.

6LetPK

k=1Ak}Bk}Ckbe the rank-K CP approximation ofG. Stack these as the columns ofA;B;C2RKK, and set [U;;V] = [U?A;?B;V?C].

7For eachk= 1;:::;K, rescale the initializers to satisfy constraints onVand.Lemma 4.1SupposeX=PK

k=1Uk}k}Vkand is completely observed. Replace? in Algorithm 2 by the topKeigenvectors ofW=1I P K k=1Y(k)>Y(k). Then, the outputs

U;;Vof Algorithm 2 satisfy thatX= (1 +)PK

k=1Uk}k}Vk.

Proof of Lemma 4.1 is given in Appendix B.2.

4.2

Auto-selection of tuning parameters

Selection of regularizers1and2:One way to choose the tuning parameters1k and2kis to use cross-validation. However, this can be computationally expensive evenquotesdbs_dbs28.pdfusesText_34

[PDF] Répertoire des cours - wwwedugovonca - Ontarioca

[PDF] code de deontologie belge francophone des assistants sociaux ufas

[PDF] code de deontologie belge francophone des assistants sociaux ufas

[PDF] code de deontologie as cpas - Comité de Vigilance en Travail Social

[PDF] code de déontologie - Ministère de l 'Éducation et de l 'Enseignement

[PDF] code de deontologie medicale algerien - ATDS

[PDF] Code de déontologie des sages-femmes - Ordre des sages-femmes

[PDF] code de deontologie des praticiens de l 'art infirmier belge

[PDF] CODE DE DéONTOLOGIE MéDICALE - Conseil National de l Ordre

[PDF] CODE DE DéONTOLOGIE MéDICALE - Conseil National de l Ordre

[PDF] Code de déontologie des médecins

[PDF] Code de déontologie des médecins

[PDF] code de deontologie medicale tunisien - ATDS

[PDF] Code de droit économique - WIPO

[PDF] Code des droits et des procédure fiscaux - ATB Entreprise