Modeling treatment events in disease progression

05/26/2019 ∙ by Guanyang Wang, et al. ∙ Stanford University 0

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.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 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 for each patient, given sparse observations . Due to the high-dimensional nature of longitudinal data, existing results usually restrict solutions to subspace of functions and utilize similarities between patients via enforcing low-rank structures. One popular approach is the mixed effect models, including Gaussian process approaches verbeke1997linear ; zeger1988models and functional principal components james2000principal . 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 to pose the problem of disease progression estimation as an optimization problem. Kidzinski and Hastie. kidzinski2018longitudinal proposed a framework which formulates the problem as a matrix completion problem and solves 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 kidzinski2018longitudinal , 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 be modeled differently.

In this work, we propose a model for tracking disease progression that includes the effects of treatments. We introduce the Coordinatewise-Soft-Impute (CSI) algorithm for fitting the model and investigate its theoretical and practical properties. The contribution of our work is threefold: First, we propose a model and an algorithm CSI, to estimate the progression of disease which incorporates the effect of treatment events. The framework is flexible, distribution-free, simple to implement 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 on both simulated data and a dataset of Gillette Children’s Hospital with patients diagnosed with Cerebral Palsy, and demonstrate the superior performances of CSI.

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 4. Finally, in Section 5 and 6 we provides empirical results of CSI on the simulated and the real datesets respectively. We discuss some future directions in Section 7.

2 Problem statement and related work

Let be the trajectory of our objective variable, such as the size of tumor, over fixed time range , and be the number of patients. For each patient , we measure its trajectory at irregularly time points and denote the results as . We are primarily interested in estimating the disease progression trajectories of all patients, based on observation data .

To fit a continuous curve based on discrete observations, we restrict our estimations to a finite-dimensional space of functions. Let be a fixed basis of (e.g. splines, Fourier basis) and be first components of it. The problem of estimating can then be reduced to the problem of estimating the coefficients such that is close to at time .

When the number of observations per patient is less than or equal to the number of basis functions , we can perfectly fit any curve without error, leading to overfitting. Moreover, this direct approach ignores the similarities between curves. Below we describe two main lines of research improving on this, the mixed-effect model and the matrix completion model.

2.1 Linear mixed-effect model

In mixed-effect models, every trajectory is assumed to be composed of two parts: the fixed effect for some that remains the same among all patients and a random effect that differs for each . In its simplest form, we assume

where is the covariance matrix,

is the standard deviation and

, are functions and evaluated at the times , respectively. Estimations of model parameters

can be made via expectation maximization (EM) algorithm

laird1982random . Individual coefficients can be estimated using the best unbiased linear predictor (BLUP) henderson1975best .

In linear mixed-effect model, each trajectory is estimated with degrees of freedom, which can still be too complex when observations are sparse. One typical solution is to assume a low-rank structure of the covariance matrix by introducing a contraction mapping from the functional basis to a low-dimensional latent space. More specifically, one may rewrite the LMM model as

where is a matrix with and 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 james2000principal ; lawrence2004gaussian ; schulam2016disease . In the later sections, we will compare our algorithm with one specific implementation named functional-Principle-Component-Analysis (fPCA) james2000principal , which uses EM algorithm for estimating model parameters and latent variables .

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 to 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 et al. in kidzinski2018longitudinal 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 into equi-distributed points with and let be the projection of the -truncated basis onto grid . The observation matrix is constructed from the data by rounding the time of every observation to the nearest time grid and regarding all other entries as missing values. Due to sparsity, we assume that no two observation ’s are mapped to the same entry of .

Let denote the set of all observed entries of . For any matrix , let be the projection of onto , i.e.  where for and otherwise. Similarly, we define to be the projection on the complement of . Under this setting, the trajectory prediction problem is reduced to the problem of fitting a matrix such that on observed indices .

The direct way of estimating is to solve the optimization problem


where is the Fröbenius norm. Again, if is larger than the number of observations for some subject we will overfit. To avoid this problem we need some additional constraints on . A typical approach in the matrix completion community is to introduce a nuclear norm penalty—a relaxed version of the rank penalty while preserving convexity rennie2005fast ; candes2009exact . The optimization problem with the nuclear norm penalty takes form


where is the regularization parameter, is the Fröbenius norm, and

is the nuclear norm, i.e. the sum of singular values. In

kidzinski2018longitudinal , a Soft-Longitudinal-Impute (SLI) algorithm is proposed to solve (2.2) efficiently. For completeness of the paper, we include the SLI algorithm in Appendix A, while noting that it is also a special case of our algorithm 1 defined in the next section with fixed to be .

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 is impossible to build a 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.1), we build our model based on matrix completion framework of Section 2.2.

More specifically, let be the time of treatment of the ’th patient, rounded to the closest ( if no treatment is performed). We encode the treatment information as a zero-one matrix , where if and only , i.e. patient has already taken the treatment by time . Each row of takes the form of . Let denote the average additive effect of treatment among all patients. In practice, we have access to the sparse observation matrix and surgery matrix and aim to estimate the treatment effect and individual coefficient matrix based on and the fixed basis matrix such that .

Again, to avoid overfitting and exploit the similarities between individuals, we add a penalty term on the nuclear norm of . The optimization problem is thus expressed as:


for some .

3.1 Coordinatewise-Soft-Impute (CSI) algorithm

Though the optimization problem (3.1) above does not admit an explicit analytical solution, it is not hard to solve for one of or given the other one. For fixed , the problem reduces to the optimization problem (2.2) with and can be solved iteratively by the SLI algorithm kidzinski2018longitudinal , which we will also specify later in Algorithm 1. For fixed , we have


where is the set of non-zero indices of . Optimization problem (3.2) can be solved by taking derivative with respect to directly, which yields


The clean formulation of (3.3) motivates us to the following Coordinatewise-Soft-Impute (CSI) algorithm (Algorithm 1): At each iteration, CSI updates from via soft singular value thresholding and then updates from via (3.3), finally it replaces the missing values of based . In the definition, we define operator as for any matrix , , where is the SVD of and is derived from the diagonal matrix . Note that if we set throughout the updates, then we get back to our base model SLI without treatment effect.

  1. Initialize

    all-zero matrix,


  2. Repeat:

    1. Compute ;

    2. Compute ;

    3. If , exit;

    4. Assign , .

  3. Output , .

Algorithm 1 Coordinatewise-Soft-Impute

4 Convergence Analysis

In this section we study the convergence properties of Algorithm 1. Fix the regularization parameter , let be the value of in the ’th iteration of the algorithm, the exact definition of which is provided below in (4.4). We prove that Algorithm 1

reduces the loss function at each iteration and eventually converges to the global minimizer.

Theorem 1.

The sequence converges to a limit point which solves the optimization problem:

Moreover, satisfies that


The proof of Theorem 1 relies on five technique Lemmas stated below. The detailed proofs of the lemmas and the proof to Theorem 1 are provided in Appendix B. The first two lemmas are on properties of the nuclear norm shrinkage operator defined in Section 3.1.

Lemma 1.

Let be an matrix and is an orthogonal matrix of rank . The solution to the optimization problem is given by where is defined in Section 3.1.

Lemma 2.

Operator satisfies the following inequality for any two matrices , with matching dimensions:



Lemma 1 shows that in the -th step of Algorithm 1, is the minimizer for function . The next lemma proves the sequence of loss functions is monotonically decreasing at each iteration.

Lemma 3.

For every fixed , the ’th step of the algorithm is given by


Then with any starting point , the sequence satisfies

The next lemma proves that differences and both converge to .

Lemma 4.

For any positive integer , we have Moreover,

Finally we show that if the sequence , it has to converge to a solution of (4.1).

Lemma 5.

Any limit point of sequences satisfies (4.1).

5 Simulation study

In this section we illustrate properties of our Coordinatewise-Soft-Impute (CSI) algorithm via simulation study. The simulated data are generated from a mixed-effect model with low-rank covariance structure on :

for which the specific construction is deferred to Appendix C. Below we discuss the evaluation methods as well as the results from simulation study.

5.1 Methods

We compare the Coordinatewise-Soft-Impute (CSI) algorithm specified in Algorithm 1 with the vanilla algorithm SLI (corresponding to in our notation) defined in kidzinski2018longitudinal and the fPCA algorithm defined in james2000principal based on mixed-effect model. We train all three algorithms on the same set of basis functions and choose the tuning parameters (for CSI and SLI) and (for fPCA) using a 5-fold cross-validation. Each model is then re-trained using the whole training set and tested on a held-out test set consisting 10% of all data.

The performance is evaluated in two aspects. First, for different combinations of the treatment effect and observation density , we train each of the three algorithms on the simulated data set, and compute the relative squared error between the ground truth and estimation ., i.e., . Meanwhile, for different algorithms applied to the same data set, we compare the mean square error between observation and estimation over test set , namely,


We train our algorithms with all combinations of treatment effect , observation rate , and thresholding parameter (for CSI or SLI) or rank (for fPCA). For each fixed combination of parameters, we implemented each algorithm times and average the test error.

5.2 Results

The results are presented in Table 1 and Figure 1. From Table 1 and the left plot of Figure 1, we have the following findings:

  1. CSI achieves better performance than SLI and fPCA, regardless of the treatment effect and observation rate . Meanwhile SLI performs better than fPCA.

  2. All three methods give comparable errors for smaller values of . In particular, our introduction of treatment effect does not over-fit the model in the case of .

  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 , the of CSI decreases from of SLI and of fPCA at to of SLI and of fPCA at .

  4. All three algorithms suffer a higher with smaller observation rate . The biggest decay comes from SLI with an average 118% increase in test error from to . 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 to estimate , we plot the relative squared error of using CSI with different observation rate in the right plot of Figure 1. As shown in Figure 1, regardless of the choice of observation rate and treatment effect , is always smaller than and most of the estimations achieves error less than . Therefore we could conclude that, even for sparse matrix , the CSI algorithm could still give very accurate estimate of the treatment effect .

Figure 1: Left: Comparisons between fPCA, SLI and CSI in estimating with different observation rates. Lines with colors red, green and blue correspond to fPCA, SLI and CSI respectively. Dotted, dashed and straight lines correspond to observation rate , and respectively. Right: Relationship between relative squared error of and treatment effect using CSI with different observation rate. Lines with colors red, green and blue correspond to observation rate , and respectively.
Observation rate
Treatment effect
Table 1: Comparisons between fPCA, SLI and CSI under different values of and .

6 Data Study

In this section, we apply our methods to real dataset on the progression of 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. Orthopaedic surgery plays a major role in minimizing gait impairments related to CP mcginley2012single . 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 a 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 set 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 schwartz2008gait , 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 as GDI values. Due to insufficiency of data, we model surgeries of different types and multiple surgeries as a single additive effect on GDI measurements following the methodology from Section 3. We test the same three methods CSI, SLI and fPCA 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 dimensions evaluated at a grid of points, with regularization parameters for CSI and SLI and rank constraints for fPCA. To ensure sufficient observations for training, we cross validate and test our models on patients with at least 4 visits and use the rest of the data as a 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/validation/test set for 40 times and compare the mean and sd of test-errors. The results are presented in Table 3 and Figure 3. Compared with the null model pMean (Column 2 of Table 3), 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 of vanilla model SLI by 7%, it also provide a stable estimation with the smallest sd across multiple selections of test sets.

mean scaled mean sd
CSI 74.28 0.62 8.90
SLI 79.92 0.66 9.22
fPCA 127.73 1.06 13.54
rMean 87.26 0.73 8.96
pMean 119.80 1.00 12.84
Figure 2: Test error on GDI dataset
Figure 3: Box plot for test errors

We take a closer look at the low-rank decomposition of disease progression curves provided by algorithms. Fix one run of algorithm CSI with

, there are 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 

3(a). The first PC recovers the general trend that gait disorder develops 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 SLI and fPCA as well.

An example of predicted curve from patient ID 5416 is illustrated in Figure 3(b) , where the blue curve represents the prediction without estimated treatment effect , green curve the final prediction and red dots actual observations. It can be seen that the additive treatment effect helps to model the sharp difference between the exam before exam (first observation) and later exams.

(a) Top 3 PCs from CSI algorithm
(b) Predicted curve of patient ID 5416
Figure 4: Low-rank decomposition of disease progression curves

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 a 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 to mention that, although we focus on 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 to more complicated settings. In our model, treatments have been characterized as the binary matrix with a single parameter . In practice, each individual may take different types 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 . 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 procedure to answer the proposed question.

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.


  • [1] Jian-Feng Cai, Emmanuel J Candès, and Zuowei Shen. A singular value thresholding algorithm for matrix completion. SIAM Journal on Optimization, 20(4), 2010.
  • [2] Emmanuel J Candès and Benjamin Recht. Exact matrix completion via convex optimization. Foundations of Computational mathematics, 9(6):717, 2009.
  • [3] Charles R Henderson.

    Best linear unbiased estimation and prediction under a selection model.

    Biometrics, pages 423–447, 1975.
  • [4] Gareth M James, Trevor J Hastie, and Catherine A Sugar. Principal component models for sparse functional data. Biometrika, pages 587–602, 2000.
  • [5] Łukasz Kidziński and Trevor Hastie. Longitudinal data analysis using matrix completion. arXiv preprint arXiv:1809.08771, 2018.
  • [6] Nan M Laird and James H Ware. Random-effects models for longitudinal data. Biometrics, pages 963–974, 1982.
  • [7] Neil D Lawrence.

    Gaussian process latent variable models for visualisation of high dimensional data.

    In Advances in neural information processing systems, pages 329–336, 2004.
  • [8] Rahul Mazumder, Trevor Hastie, and Robert Tibshirani. Spectral regularization algorithms for learning large incomplete matrices. Journal of machine learning research, 11(Aug):2287–2322, 2010.
  • [9] 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.
  • [10] Jasson DM Rennie and Nathan Srebro. Fast maximum margin matrix factorization for collaborative prediction. In Proceedings of the 22nd international conference on Machine learning, pages 713–719. ACM, 2005.
  • [11] Peter Schulam and Raman Arora. Disease trajectory maps. In Advances in Neural Information Processing Systems, pages 4709–4717, 2016.
  • [12] Michael H Schwartz and Adam Rozumalski. The gait deviation index: a new comprehensive index of gait pathology. Gait & posture, 28(3):351–357, 2008.
  • [13] Geert Verbeke. Linear mixed models for longitudinal data. In Linear mixed models in practice, pages 63–153. Springer, 1997.
  • [14] Scott L Zeger, Kung-Yee Liang, and Paul S Albert. Models for longitudinal data: a generalized estimating equation approach. Biometrics, pages 1049–1060, 1988.

Appendix A Soft-Longitudinal-Impute (SLI) algorithm

  1. Initialize all-zero matrix.

    1. Repeat:

      1. Compute

      2. If exit

      3. Assign

  2. Output

Algorithm 2 Soft-Longitudinal-Impute

Appendix B Proofs

Proof of Lemma 1.

Note that the solution of the optimization problem


is given by (see [1] for a proof). Therefore it suffices to show the minimizer of the optimization problem (B.1) is the same as the minimizer of the following problem:

Using the fact that and , we have

On the other hand

as desired. ∎

Proof of Lemma 2.

We refer the readers to the proof in [8, Section 4, Lemma 3]. ∎

Proof of Lemma 3.

First we argue that and the first inequality immediately follows. We have

Taking derivative with respect to directly gives , as desired.

For the rest two inequalities, notice that


Here the (B.2) holds because we have

(B.3) follows from the fact that . ∎

Proof of Lemma 4.

First we analyze the behavior of ,

Meanwhile, the sequence is decreasing and lower bounded by and therefore converge to a non-negative number, yielding the differences as . Hence


as desired.

The sequence is slightly more complicated, direct calculation gives


where (B.5) follows from Lemma 2, (B.6) can be derived pairing the terms according to and .

By definition of , we have


where (B.7) follows from the Cauchy-Schwartz inequality.

Combining (B.6) with (B.7), we get

Now we are left to prove that the difference sequence converges to zero. Combining (B.4) and (B.7) it suffices to prove that . We have