Modeling treatment events in disease progression
Authors: Guanyang Wang
Yumeng Zhang
Yong Deng
Xuxin Huang
Lukasz Kidzinski
Venue:  (2020)
Keywords: disease progression
treatment events
matrix completion

Abstract: Ability to quantify and predict progression of a disease is fundamental for selecting an appropriate treatment. Many clinical metrics cannot be acquired frequently either because of their cost (e.g. MRI, gait analysis) or because they are inconvenient or harmful to a patient (e.g. biopsy, x-ray). In such scenarios, in order to estimate individual trajectories of disease progression, it is advantageous to leverage similarities between patients, i.e. the covariance of trajectories, and find a latent representation of progression. Most of existing methods for estimating trajectories do not account for events in-between observations, what dramatically decreases their adequacy for clinical practice. In this study, we develop a machine learning framework named Coordinatewise-Soft-Impute (CSI) for analyzing disease progression from sparse observations in the presence of confounding events. CSI is guaranteed to converge to the global minimum of the corresponding optimization problem. Experimental results also demonstrates the effectiveness of CSI using both simulated and real dataset.

Under review as a conference paper at ICLR 2020
MODELING TREATMENT EVENTS IN DISEASE PRO- GRESSION
Anonymous authors Paper under double-bling review
ABSTRACT
Ability to quantify and predict progression of disease is fundamental for selecting an appropriate treatment. Many clinical metrics cannot be acquired frequently either because of their cost (e.g. MRI, gait analysis) or because they are incon- venient or harmful to a patient (e.g. biopsy, x-ray). In such scenarios, in order estimate individual trajectories of disease progression, it is advantageous leverage similarities between patients, i.e. the covariance of trajectories, and find a latent representation of progression. However, our knowledge, most of ex- isting methods for estimating trajectories do not explicitly account for treatment events such as surgery in-between observations, which dramatically decreases their adequacy for clinical practice. In this study, we develop first machine learning framework named Coordinatewise-Sofff fiptte (CSI) for analyzing discase pro- gression from sparse observations in the presence of confounding events. CSI is easy to implement and is guaranteed to converge to the global minimum the corresponding optimization problem. Experimental results also demonstrate the effectiveness of CSI using both simulated and real datasets.
INTRODUCTION
The course of disease progression in individual patients is one of the biggest uncertainties in medical practice. In an ideal world, accurate, continuous assessment of a patient's condition helps with prevention and treatment. However, many medical tests are either harmful or inconvenient to perform frequently, and practitioners have to infer the development of disease from sparse, noisy observations.
In its simplest form, the problem of modeling disease progressions is to fit the curve of ( E [fmin: max]] for each patient, given sparse observations := (7(t))...... (tx))). )) Due to the high- subspace of
dimensional nature of longitudinal data, existing results usually restrict solutions functions and utilize similarities between patients via enforcing low-rank structures. One popular approach is the mixed effect models, including Gaussian process approaches (Verbeke 19977 Zeger et al. 1988) and functional principal components (James et al. 2000). While generative models are commonly used and have nice theoretical properties, their result could be sensitive to the underlying distributional assumptions of observed data and hard to adapt to different applications. Another line of research is pose the problem of disease progression estimation as an optimization problem. Kidzinski and Hastie. Kidziński & Hastic (2018) proposed a framework which formulates the problem as a matrix completion problem and solve it using matrix factorization techniques. This method is distribution-free and flexible to possible extensions.
Meanwhile, both types of solutions model the natural progression of disease using observations of the targeted variables only. They fail to incorporate the existence and effect of human interference: medications, therapies, surgeries, etc. Two patients with similar symptoms initially may have different futures if they choose different treatments. Without that information, predictions can be way-off.
To the best of our knowledge, existing literature talks little about modeling treatment effect on disease progression. In Kidziński & Hastie (2018), authors use concurrent observations of auxillary variables (e.g. oxygen consumption to motor functions) to help estimate the target one, under the assumption that both variables reflect the intrinsic latent feature of the disease and are thus correlated. Treatments of various types, however, rely on human decisions and to some extent, an exogenous variable to the development of disease. Thus they need to modeled differently.

---

Under review as a conference paper at ICLR 2020
In this work, we propose model for tracking disease progression that includes the effects of treatments. We introduce the Coordinatewise- BaarinatewieeSSff-lmppeee investigate its theoretical and practical properties. The contribution of our work is threefold: First, we propose model and an algorithm CSI, to estimate the progression of disease which incorporates the effect of treatment events. The framework is flexible, distribution-ffed and generalizable. Second, we prove that CSI converges to the global solution regardless of the initialization. Third, we compare the performance of CSI with various other existing methods both simulated data and dataset of Gillette Children's is Hospital with patients diagnosed with Cerebral Palsy, and demonstrate the superior performances of CSL.
(CSI) algorithm for fitting the model and
simple to implement
The rest of the paper is organized as follows. In Section 2 we state the problem and review existing methods. Next, in Section 3 we describe the model and the algorithm. Theoretic properties of the algorithm are derived in Section Finally, in Section and we provides empirical results of CSI on the simulated and the real datesets respectively. We discuss some future directions Section
PROBLEM STATEMENT AND RELATED WORK
Let v(t) be the trajectory of our objective variable, such as the size of tumor, over fixed time range t E [imi m.a]] and N be the number of patients. For each patient 1 < i < N, we measure its trajectory (t) at n, irregularly time points t; (4,,,, .... ti,,,' and denote the results as y: = [Vs,1, Vi... = [ys (f,,1), ys (I'.,,))'. We are primarily interested in estimating the disease progression trajectories {t)] of all N patients, based on observation data ((t...)))
To fit continuous curve based discrete observations, we restrict estimations to finite- dimensional space of functions. Let {b,, E N} be a fixed basis of L2 max)) (c.g. splines, Fourier basis) and b {b; S i < K} be first K dimensions of it. The problem of estimating (t) can then be reduced to the problem of estimating the coefficients w, = [ W1 ". Wi,k such that w(b(t) is close y/ (t) at time / E t.
Though intuitive, the above method has two main drawbacks. First, when the number of observations per patient is less than equal to the number of basis functions K, can perfectly fit any curve without error, leading to overfitting. Moreover, this direct approach ignores the similarities between curves. Different patients may share similar trend of the trajectories which could potentially imporve the prediction. Below we describe two main lines of research improving on this, the mixed-effect model and the matrix completion model.
2.1 LINEAR MIFED-EFFECT MODEL
In mixed-effect models, every trajectory (f) is assumed to be composed of two parts: the fixed effect j(t) m' mt) for some m E T that remains the same among patients and a random effect w; E R that differs for each i E [1,...,N): In its simplest form, we assume W. (O 2) and y:/w; ~ N(u + o,
where E is the K X K covariance matrix, a is the standard deviation and f4 [u(tt,1,, (((,,)) (.t,,)))...((,,m,... (ttt,,;))'' )l' Jl'," B; = [b(tt,1). b(tr,2.... ".tt:)))....bt.m... )I' are functions p(t) and b(t) eval- uated at the times t;, respectively. Estimations of model parameters p, E can be made via expectation maximization (EM) algorithm (Laird & Ware 1982). Individual coefficients w; can be estimated using the best unbiased linear predictor (BLUP) (Henderson] 1975).
In linear mixed-effect model, each trajectory is estimated with |w: K degrees of freedom, which can still be too complex when observations are sparse. One typical solution is to assume low-rank structure of the covariance matrix Σ by introducing a contraction mapping A from the functional basis to a oiw-dimensional latent space. More specifically, one may rewrite the LMM model as
yı/W; N/4 + B,Aw;, o?
where A is a K x q matrix with < K and w, E IR9 is the new, shorter random effect to be estimated. Methods based on low-rank approximations are widely adopted and applied in practice and different algorithms on fitting the model have been proposed (James et al.) 2000 Lawrence 2004] Schulam & Arora 2016). In the later sections, will compare our algorithm with one specific implementation named functional- Aimetional-Principle-- Principle-Compoeent uses EM algorithm for estimating model parameters and latent variables w;,
(IPCA) (James et al., 2000), which
2

---

Under review as a conference paper at ICLR 2020
2.2 MATRIX COMPLETION MODEL
While the probabilistic approach of mixed-effect models offers many theoretical advantages including convergence rates and inference testing, it is often sensitive the assumptions on distributions, some of which are hard to verify in practice. To avoid the potential bias of distributional assumptions in mixed-effect models, Kidzinski and Hastie (Kidziński & Hastie 2018) formulate the problem as a sparse matrix completion problem. We will review this approach in the current section.
To reduce the continuous-time trajectories into matrices, we discretize the time range [fmmin f into T equi-distributee points G [TI...., Tr] with TI = fmin: TT = fmax and let B [b(r)), ((T))..... b( E R K be the projection of the K-truncated basis b onto grid G. The N x K observation matrix Y is constructed from the data ((t; y+)] by rounding the time li.j of every observation y:((,,) to the nearest time grid and regarding all other entries as missing values. Due to sparsity, we assume that no two observation (ttt)))'s
mapped to the same entry of Y.
Let f2 denote the set of all observed entries of Y. For any matrix A. let (r((A) be the projection of A onto St, i.e. (AA = M where Mij = Ai for (i,j) E Z and Maj = 0 otherwise. Similarly, we define (2() A (A(() to be the projection on the complement of Q. Under this setting, the trajectory prediction problem is reduced to the problem of fitting a N x K matrix W such that WB' R Y on observed indices ?
The direct way of estimating W is to solve the optimization problem
1 2 |1 1B
arg min W
(2.1)
where -|-
|| the Fröbenius norm, Again, if K is larger than the number of observations for some subject we will overfit. To avoid this problem we need some additional constraints on W. A typical approach in the matrix completion community is to introduce nuclear norm penalty- a relaxed version of the rank penalty while preserving convexity (Rennie & Srebro 2005 Candès & Recht 2009). The optimization problem with the nuclear norm penalty takes form
arg min min W
3 P1 WB) + \||
(2.2)
where À > 0 is the regularization parameter, ||1|| p the Fröbenius norm, and || I|. is the nuclear norm, i.e. the sum of singular values. In Kidzinski & Hastie (2018), Soft-Longitudinll-nnn (SLI) algorithm is proposed to solve (2.2) efficiently. We refer the readers to Kidziński & Hastie
(2018) for detailed description of SLI while noting that it is also special case of our algorithm defined in the next section with fixed to be O.
3 MODELING TREATMENT IN DISEASE PROGRESSION
In this section, we introduce our model on effect of treatments in disease progression.
A wide variety of treatments with different effects and durations exist in medical practice and it impossible to build single model to encompass them all. In this study we take the simplified approach and regard treatment, with the example of one-time surgery in mind, as a non-recurring event with an additive effect on the targeted variable afterward. Due to the flexibility of formulation of optimization problem (2.I), we build our model based on matrix completion framework of Section 2.2
More specifically,, let s(i) G be the time of treatment of the i'th patient, rounded to the closest Th G (s(1) = co if no treatment is performed). We encode the treatment information as N x T zero-one matrix Ts. where (Ts) = if and only Tj 2 si)), i.e. patient has already taken the treatment by time Ty. Each row of Is takes the form of (0.... ... 0, 1, 1). Let μ denote the average additive effect of treatment among all patients. In practice, we have access to the sparse observation matrix Y and surgery matrix Ig and aim to estimate the treatment effect μ and individual coefficient matrix W based on Y, Is and the fixed basis matrix B such that WB' + pIs = Y.
Again, to avoid overfitting and exploit the similarities between individuals, we add a penalty term on the nuclear norm of W. The optimization problem is thus expressed as:
arg min 1 / P1 WB' - -Is)ll + X||W]|+
(3.1)
M.W 2
3

---

Under review as a conference paper at ICLR 2020
for some À 0.
3.1 COORDINATEWISE-SOFFE SOFTIISPSSS (CST) ALGORITHM
Though the optimization problem (3..) above does not admit an explicit analytical solution, is not hard to solve for of μ W given the other one. For fixed J4, the problem reduces to the optimization problem (.2) with Y Y pls and can be solved iteratively by the SLI
algorithm Kidziñski & Hastie 2018), which we will also specify later in Algorithm For fixed W, have
arg min |PPa(Y - WB' - +
1 2
1 arg min _ 1
- = arg min
1
Σ
(YY - k, - -4),
(3.2)
2
P
where Ggg is the set of non-zero indices of Is. Optimization problem 3.2 can be solved by taking derivative with respect to directly, which yields
N Eo.nenna (Y - WB
(3.3)
[an
The clean formulation of (.3 motivates to the following
wwwpatttwwiiSSSt---ppeeeee (CSI) algorithm (Algorithm D:: each iteration, CSI updates Now from (Wous: Hoia via soft singular value thresholding and then updates Mnew from (Wnew, (Mola) via (333. finally it replaces the missing values of Y based (Wnew: /.ew). In the definition, we define operator Sx as for any matrix X, (x))) = UDVV, where X = DDV is the SVD of X and D> = diag((max - ()) is derived from the diagonal matrix D diag((d;)) 1). Note that if we set = 0 throughout the updates, then we get back to our base model SLI without treatment effect.
Algorithm 1: SAAILIAATEEIII-SSS--MMPPPES
1. Initialize Wald f all-zero matrix, Hold + O.
2. Repeat:
(a) Compute W new SX(PPo(Y Molols) Pt (WoaB"))B
Eens)eeenen (Y- Wan B'),,
(b) Compute Pnew +
(c) max { (Mezn m | Woull? uāu
2 < exit:
(d) Assign Wold & W Wnews Hold + Hnew-
3. Output WA Wnew, jex Unew-
4 CONVERGENCE ANALYSIS
In this section we study the convergence properties of Algorithm Fix the regularization parameter
> 0, let (1L) (k) (k) be the value of (μ,, W) in the k'th iteration of the algorithm, the exact definition of which is provided below in (4.4). We prove that AlgorithmTreeaae Treduces the loss function at
each iteration and eventually converges to the global minimizer.
Theorem 1. The sequence (K) ((()) converges limit point (ix, Wx) which solves the opti- mization problem:
(Ax, WW) arg min n.W
Pa(Y WB' -1s)) + +
Moreover, (fx, Wx) satisfies that
Axls) + P P Ax E Ensennng -Y- -WWBBka
WA = SX((Po
(4.1)
4

---

Under review as a conference paper at ICLR 2020
The proof of Theorem I relies five technique Lemmas stated below. The detailed proofs of the lemmas and the proof to Theorem are provided in Appendix [A] The first two lemmas are on properties of the nuclear norm shrinkage operator Sx defined in Section 3.1
Lemma 1. Let W be N x K matrix and B is orthogonal T x K matrix of rank K. The solution to the optimization problem minw 2 WB + W| is given by W Sx B) where Sx B) is defined Section 3.7
Lemma 2. Operator Sx(+) satisfies the following inequality for any two matrices W,, W2 with
matching dimensions:
(W)) Sx(W2) & [ W2
Define
P. = 2 P1 WB' - pls) A/[W,
(4.2)
1 11PPa pTs) + (WBB WB + AI[W]|.
Qx W. p)
(4.3)
Lemma shows that in the k-th step of Algorithm Qx(- (x(-1W((-1), ()))),) The next lemma proves the sequence of loss functions (W MA is monotonically decreasing at each iteration.
(( is the minimizer for function
Lemma 3. For every fixed 2 O, the k k'th step of the algorithm (ux) ((())) is given by
Easeanas (Y xx()) B')ss (4.4)
(( (k) arg min Q QW
Px 11
(K)
Then with any starting point (1) ())) the sequence {
(k) (()
} satisfies
(* (44)) S ( Px
(k- (( s (W (((-))) 1x (k- 1)) s fa (W Mx (1-1))).
The next lemma proves that differences (μt - 1k-1)² and [|w 1/ || both converge to O.
Lemma 4. For any positive integer we have [|w-(kt1) w/(h) || s | ( 12. Moreover,
(k+1)
(K)
O, W w((h) +0 as
11 ux
+ oc.
Finally we show that if the sequence (1() (x) it has to converge to solution of 4.1
}a,
Lemma 5. Any limit point (Axx Wx) of sequences {(p) ((()))> } satisfies 4.1
SIMULATION STUDY
In this section we illustrate properties of our Coordinatewise-Soff Impuut (CSI) algorithm via simulation study. The simulated data are generated from mixed-effect model with low-rank
covariance structure on
Y WB pls +
for which the specific construction is deferred Appendix B Below discuss the evaluation methods as well as the results from simulation study.
5.1 METHODS
We compare the Coordinatewise- Soft Soft-Immute (CSI) algorithm specified in Algorithm with the vanilla algorithm SLI (corresponding to u = 0 in our notation) defined in Kidziński & Hastie (2018) and the fPCA algorithm defined in James et al. (2000) based on mixed-effect model. We train all three algorithms the same set of basis functions and choose the tuning parameters À (for CSI and
5

---

Under review as a conference paper at ICLR 2020
SLI) and R (for fPCA) using 5-fold cross-validation. Each model is then re-trained using the whole training set and tested on a held-out test set Qusa consisting 10% of all data.
The performance is evaluated in two aspects. First, for different combinations of the treatment effect and observation density p. we train each of the three algorithms on the simulated data set, and compute the relative squared error between the ground truth f and estimation fa., i.e., RSE(A) ( n1)² / Meanwhile, for different algorithms applied to the same data set, we compare the mean square error
between observation Y and estimation Y over test set Sreat. namely,
MSE(V)
Σ
(Vis - Viz)² = [Ress Pa (Y) - Pa (
(5.1)
(4.5) EDBon
We train our algorithms with all combinations of treatment effect E {0, 00.00.00.4.-... 5). obser- vation rate p E {0.1. 0.5). and thresholding parameter À E (0, 4]
(for CSI or SLI) or rank R {2.3, ... 6} (for fPCA). For each fixed combination of parameters, implemented each
algorithm 10 times and average the test error.
5.2 RESULTS
The results are presented in Table I and Figure From Table and the left plot of Figure 1 the following findings:
we have
1. CSI achieves better performance than SLI and fPCA, regardless of the treatment effect j and observation rate p. Meanwhile SLI performs better than fPCA.
2. All three methods give comparable errors for smaller values of μ. In particular, our introduc- tion of treatment effect μ does not over-fit the model in the case of J O.
3. As the treatment effect increases, the performance of CSI remains the same whereas the performances of SLI and fPCA deteriorate rapidly. As a result, CSI outperforms SLI and fPCA by a significant margin for large values of For example, when p 0.1, the MSE(Ý) of CSI decreases from 72.3% of SLI and 59.6% of fPCA at 1 to 12.4% of SLI and 5.8% of (PCA at
4. All three algorithms suffer a higher MSE(Ý) with smaller observation rate p. The biggest decay comes from SLI with an average 118% increase in test error from p 0.5 to p= 0.1. The performances of fPCA and CSI remains comparatively stable among different observation rate with a 6% and 12% increase respectively. This implies that our algorithm is tolerant to low observation rate.
To further investigate CSI's ability estimate μ4, we plot the relative squared error of j using CSI with different observation rate in the right plot of Figure [ As shown in Figure regardless of the choice of observation rate p and treatment effect μı, RSB( always smaller than 1% and most of the estimations achieves error less than 0.1%. Therefore we could conclude that, even for sparse matrix Y. the CSI algorithm could still give very accurate estimate the treatment effect μ.
Observation rate p
0.1
0.3
0.5
Treatment effect μ
1
2
5
5
1
2
5
IPCA
0.521
2.172
5.455
0.525
2.039
5.170
0.525
2.036
5.166
MSE(Y) SLI
0.430
1.162
2.561
0.379
0.658
1,203
0.341
0.543
0.893
CSI
0.311
0.306
0.318
0.314
0.297
0.320
0.294
0.299
0.295
Table 1: Comparisons between fPCA, SLI and CSI under different values of p and p.
6
DATA STUDY
In this section, we apply our methods to real dataset on the progression motor impairment and gait pathology among children with Cerebral Palsy (CP) and evaluate the effect of orthopaedic surgeries.
Cerebral palsy is a group of permanent movement disorders that appear in early childhood. Or- thopaedic surgery plays major role in minimizing gait impairments related to CP (MeGinley et al.

---

Under review as a conference paper at ICLR 2020
0.0020
Mothods
Incan
IPCAD-X
0.005..
SPCAD.S
Oservatton rate
56.011
0.0010
2-
C310.3
000s
C919.5
0.0000
Treatment offect 14
Treatment affectt
Figure 1: Left: Comparisons between fPCA, SLI and CSI in estimating Y with different obser- vation rates. Lines with colors red, green and blue correspond fPCA, SLI and CSI respectively. Dotted, dashed and straight lines correspond to observation rate P 0.1, 0.3 and 0.5 respectively. Right: Relationship between relative squared error of ji and treatment effect ju using CSI with different observation rate. Lines with colors red, green and blue correspond to observation rate 0.1, 0.3 and 0.5 respectively.
2012). However, it could be hard to correctly evaluate the outcome of a surgery. For example, the seemingly positive outcome of a surgery may actually due to the natural improvement during puberty, Our objective is to single out the effect of surgeries from the natural progression of disease and use that extra piece of information for better predictions.
6.1 DATA AND METHOD
We analyze data set of Gillette Children's Hospital patients, visiting the clinic between 1994 and 2014, age ranging between 4 and 19 years, mostly diagnosed with Cerebral Palsy. The data contains 84 visits of 36 patients without gait disorders and 6066 visits of 2898 patients with gait pathologies. Gait Deviation Index (GDI), one of the most commonly adopted metrics for gait functionalities (Schwartz & Rozumalski 2008 was measured and recorded at each clinic visit along with other data such as birthday, subtype of CP, date and type of previous surgery and other medical results.
Our main objective is to model individual disease progression quantified GDI values. Due to insufficiency of data, we model surgeries of different types and multiple surgeries as single additive effect on GDI measurements following the methodology from Section 3 We test the same three methods CSI, SLI and IPCA as in Section 5 and compare them to two benchmarks the population mean of all patients (pMean) and the average GDI from previous visits of the same patient (rMean).
All three algorithms was trained on the spline basis of K 9 dimensions evaluated at a grid 7 51 points, with regularization parameters À E {20, 25, .... 40) for CSI and SLI and rank constraints r (2...., 6} for iPCA. To ensure sufficient observations for training, we cross validate and test models patients with at least 4 visits and use the rest of the data common training set. The effective size of 2-fold validation sets and test set are 5% each. We compare the result of each method/combination of parameters using the mean square error of GDI estimations on held-out entries as defined in (5.1)
6.2 RESULTS
We run all five methods on the same training/validationnesst set for 40 times and compare the mean and sd of test-errors. The results are presented in Table Zan Figure Compared with the null model MMen (Column of Table fPCA gives roughly the same order of error; CSI, SLI and rowMean provide better predictions, achieving 62%, 66% and 73% of the test errors respectively. In particular, our algorithm CSI improves the result vanilla model SLI by 7%, it also provide a stable estimation with the smallest sd across multiple selections of test sets.
We take a closer look at the low-rank decomposition of disease progression curves provided by algorithms. Fix run of algorithm CSI with X, 30, there 6 non-zero singular value vectors, which we will refer as principal components. We illustrate the top 3 PCs scaled with corresponding singular values in Figure 3a The first PC recovers the general trend that gait disorder develops
7

---

Under review as a conference paper at ICLR 2020
mean
scaled mean
sd
150-
CSI
74.28
0.62
8.90
Test Aamen
SLI
79.92
0.66
9.22
IPCA
127.73
1.06
13.54
Mean
87.26
0.73
8.96
70-
PMean
119.80
1.00
12.84
S11
PPCA Methods
pileon
Table Test error on GDI dataset
Figure 2: Box plot for test errors
76-
70-
group
eugeryy Eugey
-100
60-
10
Age
(a) Top PCs from CSI algorithm
(b) Predicted curve of patient ID 5416
Figure 3: Low-rank decomposition of disease progression curves
through age 1-10 and partially recovers during puberty. The second and third PC reflects fluctuations during different periods of child growth. By visual inspection, similar trends can be find in the top components of SL.I and fPCA as well.
An example of predicted curve from patient ID) 5416 is illustrated in Figure 3b where the blue curve represents the prediction without estimated treatment effect f = 4.33, green curve the final prediction and red dots actual observations. It can be seen that the additive treatment effect helps model the sharp difference between the exam before exam (first observation) and later exams.
7 CONCLUSION AND FUTURE WORK
In this paper, we propose a new framework in modeling the effect of treatment events in disease progression and prove corresponding algorithm CSI. To the best of our knowledge, it's the first comprehensive model that explicitly incorporates the effect of treatment events. We would also like mention that, although we focus the case of disease progression in this paper, our framework is quite general and can be used to analyze data in any disciplines with sparse observations as well as external effects.
There are several potential extensions to our current framework. Firstly, our framework could be extended more complicated settings. In our model, treatments have been characterized as the binary matrix Is with a single parameter of surgeries for one or multiple times. Secondly, the treatment effect may be correlated with the latent variables of disease type, and can be estimated together with the random effect uy Finally, our framework could be used to evaluate the true effect of a surgery. A natural question is: does surgery really help? CSI provides estimate of the surgery effect it would be interesting to design certain statistical hypothesis testing/casual inference procedure to answer the proposed question.
In practice, each individual may take different types
Though we are convinced that our work will not be the last word in estimating the disease progression, we hope our idea is useful for further research and we hope the readers could help to take it further.
REFERENCES
Jian-Feng Cai, Emmanuel Candès, and Zuowei Shen. A singular value thresholding algorithm for matrix completion. SIAM Journal Optimization, 20(4), 2010.
8

---

Under review as a conference paper at ICLR 2020
Emmanuel J Candès and Benjamin Recht. Exact matrix completion via convex optimization. Foun- dations of Computational mathematics, 966):11,, 2009.
Charles R Henderson. Best linear unbiased estimation and prediction under a selection model. Biometrics, pp. 423-447, 1975.
Gareth James, Trevor Hastie, and Catherine Sugar. Principal component models for sparse functional data. Biometrika, pp. 587-602, 2000.
Lukasz Kidziński and Trevor Hastie. Longitudinal data analysis using matrix completion. arXiv preprint arXiv:1809.0877) 2018.
Nan M Laird and James H Ware. Random-effects models for longitudinal data. Biometrics, pp. 963-974, 1982.
Neil D Lawrence. Gaussian process latent variable models for visualisation of high dimensional data. In Advances in neural information processing systems, pp. 329-336, 2004.
Rahul Mazumder, Trevor Hastie, and Robert Tibshirani. Spectral regularization algorithms for learning large incomplete matrices. Journal of machine learning research, 1((Aug)2287-2322 2010.
Jennifer L McGinley, Fiona Dobson, Rekha Ganeshalingam, Benjamin J Shore, Erich Rutz, and H Kerr Graham. Single-event multilevel surgery for children with cerebral palsy: a systematic review. Developmental Medicine & Child Neurology, 54(2):117-128, 2012.
Jasson DM Rennie and Nathan Srebro. Fast maximum margin matrix factorization for collaborative prediction. In Proceedings of the 22nd international conference on Machine learning, pp- 713-719. ACM, 2005.
Peter Schulam and Raman Arora. Disease trajectory maps. In Advances in Neural Information Processing Systems, pp. 4709-4717, 2016.
Michael H Schwartz and Adam Rozumalski. The gait deviation index: comprehensive index of gait pathology. Gait & posture, 28(3):351-357, 2008
Geert Verbeke. Linear mixed models for longitudinal data. In Linear mixed models in practice, pp. 63-153. Springer, 1997.
Scott L Zeger, Kung-Yee Yee Liang, and Paul S Albert. Models for longitudinal data: generalized estimating equation approach. Biometrics, pp. 1049-1060, 1988.

---

Under review as a conference paper at ICLR 2020
A
PROOFS
Proof of Lemmal
Note that the solution of the optimization problem
min 1 Z
All + A||
(A.1)
is given by A = (x( (see Cai et al. 2010) for proof). Therefore it suffices to show the minimizer of the optimization problem (A.1) is the same as the minimizer of the following problem:
1
min W
2
Y w + I
Using the fact that ||A|12- = Tr(AA') and B'B B-- = Ir we have
W + W||W|| arg min JYB
arg min
(Tr(Y + Tr(WW") 2T((YBW')) + A|lll--
2
arg min W
(Tr(WW') - 2TT(Y + AIIII.
2
On the other hand
arg min |Y WB'B + W|ll
arg min
(Tr(YY') + Tr(WW' - 2T((BWW")) + I|Ml--
arg min
(r(WW'') 2Tr(YBW')) *||l||+
arg min
YYB - W A
2
(Y
as desired.
Proof of Lemma 2 3).
We refer the readers the proof in Mazumder et al. (2010) Section 4, Lemma
Proof of Lemma
First argue that J arg min, fx(W(()) .14) and the first inequality imme-
diately follows.
We have
arg min fa (a(W .14) = arg min Pa((Y (W(x) B' -1s)
Σ
((Y ((()) 3's 12)
arg min
(nn))nnns
Taking derivative with respect to p directly gives P = arg min, SA (W)())) .1)), as desired. For the rest two inequalities, notice that
(k) .1 (k-1))
2 Paa(Y (x( B' 1 13)1 (
1
2 PPa(Y uA (k- Is) + P
1
w BTI + (1 (A.2)
(x( (W ,11 (k-1))
S QA(W (K-1) 1-1) W (k-1) .14 (1-1))
(A.3)
11 (k- B¹
(k-1)
(k-1)
1s)
+N|lW((-1))
=A(W((E-)) .14 (k-1)).
10

---

Under review as a conference paper at ICLR 2020
Here the A.2
holds because have
マ PllPa(Y Mn (() Psi (W B¹) (1(1) 13
PIlPa MX Is /1 B¹) + Ps (W (1 B' - w B
(k-1)
2 PPa(Y - -M 15 5 B3) |/ + 1 (W (k-)) B² - B B²)
(k-1)
1 2/11 - -MX
(k-1) 11 Is ( B3) 1/.
follows from the fact that ((() arg minw Qx (WIW (k-1) :14 1)
A.3
Proof of Lemma
First we analyze the behavior of (1) J.
(k)
(A .2 -
(a x/ = - P11P - - B' - PX 1s)1l
(k-1)
1PPP( - Wx B' 1 1s)il
ssl (12 14 12-1
(B)
), fa (W creasing and lower bounded by 0 and therefore converge to non-negative number, yielding the (k)
Meanwhile, the sequence (.... fa fa .2X ( (WW .A ((
(k)
X 14 (1t) ....) is de-
)...)
differences (x( uu ( (W .u² (1))) + 0 as k - cc. Hence
(R)
(b) 12
(k-1) +
(A.4)
as desired.
The sequence (w) } is slightly more complicated, direct calculation gives
(W 1- B') ( (M) (- 11 1| = [|Sx(Pa(Y 11 (k-1) [s) (s) +
SX(Pa - -11 + (W (1-2 B'))
$ 11Pa(Y - -M² -H(t1))s) 1s) +P B)))
(A.5)
-PalY -/ - (W
(k- B'
-
= =120 (k-1) ( PP B (A.6)
where A.5 follows from Lemma A.6 can be derived pairing the 4 terms according to Ps and Pa.
(k)
By definition of ft we have
2
(k-1) 12nΩsl(u Mx Anasl
1
Σ
(W B'- -
(iJ)ERnns
(( ( B'
(A.7)
where (.7 follows from the Cauchy-Schward inequality.
Combining A.6 with (.7). we get
(k)
w (A-
1)
W
2) B'1?? = 11 - Wr (k)) (k 1) B' (x(
Now we are left to prove that the difference sequence (Wr - (k) (1-111) convverees to zero Combining
A.4 and (A77) it suffices to prove that || (H>(k-1) B' w(6 B') B¹) - 0. We have
fx( 1) .1A
[k- 2)) (x(W((--))))) ( 2) pp (k-2)) Pl (W (k 11 B' w X
2) B'
A

---

Under review as a conference paper at ICLR 2020
and the left hand side converges to because
ž fA .A2 2) Qx | ( 1112 2
(
2 faW (K-2) ( -AAW
pa ( + 0,
.H2
which completes the proof.
Proof of Lemma 5 Let (A) (ma)) w(lm)) + (ūx, WA). then Lemma gives (12 (m==1) w(mee
(px: Wx). Since we have
w(ma) Sx((Pa(( - 11 (mL
Pot (W (1 B)) E(ijennna (Y 1mm=-1) B))
1s) (
(m²)
MX
Taking limits on both sides gives us the desire result.
Proof of Theorem 7 Let (fx, Wx) be one limit point then we have:
IWa- (>(+ || = ||Sx((Pa(Y -alls) + PD (W
(A.8)
-SA((Pa -1 - ( B B)
I <snssl M4 IIP W/ ( 1'))
(A.9)
here A.8 uses Lemma 5 and A.9 uses Lemma] Meanwhile,
12 Eo.iennn (WW,
- 11 - (y(K-D))B") ( (a(k-11) 1) BB !g
12 1xx
(A.10)
Combining (A.9) and (A.10), we have
|| ( ||r & l|Ws - w(
Hence the sequence ||V/ - (x) || p monotonically decreasing and has a limit, But since there exists (mmm2) converging to Wx. the limit equals 0, which proves wi(k) W,, 1x f..
(k)
Therefore we have proved the sequence (u) ((()) always converges.
Mcanwhile, the first part of equation 4.1 and Lemma 5 in Mazumder
2010) guarantees
Ō E Ow fa (Wx, fi). By taking derivative directly we have = f. P). Therefore (Wx, A) is stationary point for (W, μ). Notice that the loss function fa 14) is a convex function with respect
to (W, p). Thus have proved that the limit point (Wx, A) minimizes the function fx p)) O
B DATA GENERATION
Let G be the grid of I' equidistributed points and let B be the basis of K spline functions evaluated on grid G. We will simulate the x observation matrix Y with three parts
Y= WB B + pls E
where W follows with uniformly distributed starting time and & represents the i.i.d. measurement error. The specific procedures is described below.
distribution with low rank structure, Is is the treatment matrix
1. Generating W given parameters x E (0, (0,1) E R, E ERK. 20
(a) Sample two K x K orthogonal matrices Vi, 1/2 via vngeular-vlte-deeoomossngg two random matrix.
(b) Sample unit length K vectors Fi, 52 via normalizing i.i.d. normal samples.
12

---

Under review as a conference paper at ICLR 2020
(c) Draw vector 8 E RN from i..d. Bernouli(x:) samples. Denote the all one vector by I.
(d) Draw N K matrices U1, U2 from i.i.d. standard normal random variables.
(e) Set
W - r. [r + U1 diagl Vil - [227
U2 diagly / V/]
where diag[s] is the diagonal matrix with diagonal elements s, "." represents coordi- natewise multiplication, and we are recycling E, i and ry7 to match the dimension.
2. Generating Is given parameter pir E (0, 1).
(a) For each 4: 1, . N, sample Tę uniformly random from {1,.... [T/pa])
(b) Set Is (1{j 2 T.sssss
3. Given parameter e Rgo, & drawn from from i.i.d. Normal(O, e²) samples.
4. Given parameter jr R, let Yo + WB pls &
5. Given parameter p (0, 1), drawn 0-1 matrix In from ii.d. Bemoulli(p) samples. Let Ω denote the set of non-zero entries of Io namely, the set observed data. Set
Y (YG)IsISNISSTT where Y.j -(o))s if (Ja)i =1
NA
otherwise
In actual simulation, fix the auxiliary parameters follows,
K = 7,7 = 51, N = 500,
K 0.33, r 1, 2,
[1,0.4, 0.005, 0.1 exp( - ..., 0.1 exp(- (app((KK + 1)], 82
0.2, 0..3 0.005, 0.1 exleex((3) ..., 0.1 exleep(-KK + 1)],
Pr 0.8, 0.5.
The remaining parameters are treatment effect μ and observation rate p. which we allow to vary across different trials.
13


Reviews:
[{"review_id": "rkgJBcaXYr", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 1, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper proposes to incorporate treatment effect in disease progression modeling for individual patients. It views treatment effect as an exogenous variable and represent that with one additional parameter in the model.\n\nThere is not enough novelty with the proposed method. It only extends Kidzinski and Hastie (2018) by adding one extra parameter that represents the average additive treatment effect. And this new model does not outperform the previous model in  Kidzinski and Hastie (2018) significantly.\n\nThe method also does not perform significantly better on real data. Given the real data  as a 1-d vector of length 6k, the method learns a 9x9 latent matrix to represent the data. And the method beats, but not significantly, the intuitive weak baseline rMean which uses the mean of the same patient’s previous visits for future prediction.\n\n\n\n\n", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I did not assess the derivations or theory.", "review_assessment:_checking_correctness_of_experiments": "I did not assess the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper."}, "extra_text": {}}, {"review_id": "S1xrtQZTKr", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 1, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper introduces Coordinatewise-Soft-Impute (CSI) algorithm that attempts to estimate the progression of a disease while taking into account the effect of a single binary treatment intervention. The authors also provide prove that CSI converges to the global solution regardless of the initialization. The paper is concluded with experimental results on synthetic  and real-world data, comparing the proposed method with contender methods functional Principal Component Analysis (fPCA) and Soft-Longitudinal-Impute (SLI).\n\nThis paper should be rejected due to the following arguments:\t\n\t- The authors claim that the “existing literature talks little about modeling treatment effect on disease progression.” However, this is not true. Here are a few works and the references within, for instance:\n\t\tXu, Y., Xu, Y., & Saria, S. (2016, December). A Bayesian nonparametric approach for estimating individualized treatment-response curves. In Machine Learning for Healthcare Conference (pp. 282-300).\n\t\tLim, B. (2018). Forecasting treatment responses over time using recurrent marginal structural networks. In Advances in Neural Information Processing Systems (pp. 7483-7493).\n\n\t- The proposed method is simple extension of that of (Kidzinski & Hastie, 2018), by adding a term that accounts for treatment and its timing. It does not present enough novelty.\n\n\t- CSI is empirically compared with fPCA and SLI, however, neither of these two methods take into account the effect of applying treatment. This comparison is unfair and therefore, cannot justify the superiority of CSI.\n\nThings to improve the paper that did not impact the score:\n\t- Citations throughout the paper had wrong formatting. Often \\citep should have been used instead of \\citet.\n\t- Page 1, last par., last sentence: Thus, they need to “be” modeled differently.\n\t- Page 2, last par., line -2: Principle >> Principal\n\t- Page 3, par. 2, line 4: the observation matrix Y is N \\times T (not N \\times K).\n\t- Page 4, line 2 of Eq. (3.2): first argument misses Y.\n\t- Code was included but no description or readme file was provided for how to run the experiments.", "extra_scores": {"experience_assessment": "I have published one or two papers in this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_experiments": "I carefully checked the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I carefully checked the derivations and theory."}, "extra_text": {}}, {"review_id": "BJlAL5MN9S", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 1, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper proposed a framework for modeling disease progression taking into account the treatment events.\n\nComments:\n1. There is no significant theoretical innovation in this paper. Specifically, the learning approach for optimizing the parameter set W in the presence of \\mu is not much different from the one proposed by (Kidzinski & Hastie, 2018) and the optimal value of \\mu is simply the solution for L2 loss function. Besides, the narrative of section 2 shares too much similarity with the context of Kidzinski & Hastie, 2018.\n\n2. The practical significance of the paper is rather limited considering that the framework can only take into account one specific treatment at a time when training a model. There are usually different treatments for a certain disease and it is hard to control enough patients to all take the same treatment. Sometimes, even the same patient could take multiple treatments at the same time or along his trajectory. \n\n3. The paper is poorly written. There are too many grammatical errors. Some sentences are hard to understand. E.g. “Treatments of various types, however, rely on human decisions and to some extent, an exogenous variable to the development of disease. Thus they need to modeled differently.” should be “Treatments of various types, however, rely on human decisions and to some extent, are exogenous variables to the development of disease. Thus they need to be modeled differently.”. “we add a penalty term on the nuclear norm of W .”, Do you mean “we add a penalty term which is the nuclear norm of W .”? Cause I don’t see any additional penalty term in the nuclear norm compared with equation 2.2. “...where (I_S)_{i,j} = 1 if and only T_i ≥ s(i)...” should be “...where (I_S)_{i,j} = 1 if and only if T_i ≥ s(i)...”. “Finally we show that if the sequence {(\\mu_{\\lambda}^{(k)}, W_{\\lambda}^{(k)} )}_k, it has to converge to a solution of (4.1).” Grammatically speaking, the if clause is not even a complete clause. ", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The proposed method has very weak novelty.
Decision: Reject