Analysing disability course in MS

This tutorial illustrates how to use the pymsprog package to analyse the evolution of disability in multiple sclerosis (MS) based on repeated assessments through time of an outcome measure (EDSS, NHPT, T25FW, SDMT; or any custom outcome measure). We’ll start by illustrating the type of input data needed, and by giving a minimal working example to introduce the main function and facilitate the reading of the document. We’ll then move on to a more detailed description of the different parameter configurations.

import pymsprog
print(pymsprog.__version__)
1.0.7

Input data

The data must be organised in a pandas DataFrame containing (at least) the following columns:

  • Subject IDs

  • Visit dates

  • Outcome values.

The visits should be listed in chronological order (if they are not, MSprog() will sort them automatically before analysing them).

For relapsing-remitting MS patients, an additional DataFrame with the dates of relapses is needed to correctly assess progression and characterise progression events as relapse-associated or relapse-independent. The dataset should contain (at least) the following columns:

  • Subject IDs

  • Visit dates.

In this tutorial, we will use toy data with artificially generated EDSS and SDMT assessments and relapse dates for a few hypothetical patients:

from pymsprog import load_toy_data
toydata_visits, toydata_relapses = load_toy_data()

print('\nVisits:')
print(toydata_visits.head())
print('\nRelapses:')
print(toydata_relapses.head())
Visits:
   id        date  visit_day  EDSS  SDMT
0   1  2021-09-23          0   5.5    54
1   1  2021-11-03         41   5.5    54
2   1  2022-01-19        118   5.5    57
3   1  2022-04-27        216   5.5    55
4   1  2022-07-12        292   6.0    57

Relapses:
   id        date  visit_day
0   2  2021-06-12        198
1   2  2022-10-25        698
2   3  2022-12-01        409
3   6  2022-12-18        426
4   7  2021-09-11        185

Dates are interpreted as calendar days (e.g., YYYY-mm-dd) by default, but can optionally be provided as number of days “from start” (starting point can be different across subjects – e.g., days from randomisation in a clinical trial): see visit_day column in the toy datasets, and SDMT example below. Disability assessment dates and relapse dates must be given in the same format.

Minimal example

The core block of pymsprog is the MSprog() function. Given data on visits and relapses in the form specified above, MSprog() detects the confirmed disability worsening (CDW) or improvement (CDI) events for each subject for the outcome of interest. If specified, CDW events can be further classified based on their timing with respect to relapses and identified as relapse-associated worsening (RAW) or progression independent of relapse activity (PIRA) [Kappos et al., 2020, Lublin et al., 2014].

By default, MSprog() detects the first CDW.

We can test it on the EDSS toy data…

from pymsprog import MSprog
summary_edss, results_edss = MSprog(toydata_visits,  # data on visits
                         subj_col='id', value_col='EDSS', date_col='date',  # specify column names
                         outcome='edss',  # specify outcome type <--- EDSS
                         relapse=toydata_relapses) # data on relapses
---
Outcome: edss
Confirmation over: 84 (-7, +730.5) days
Baseline: fixed
Relapse influence (baseline): [30, 0] days
Relapse influence (event): [0, 0] days
Relapse influence (confirmation): [30, 0] days
Events detected: firstCDW
---
Total subjects: 7
---
Subjects with CDW: 4

...or on the SDMT toy data:
summary_sdmt, results_sdmt = MSprog(toydata_visits,  # data on visits
                         subj_col='id', value_col='EDSS', date_col='visit_day',  # specify column names
                         outcome='sdmt',  # specify outcome type <--- SDMT
                         relapse=toydata_relapses,  # data on relapses
                         date_format='day')  # specify that dates are given as days
---
Outcome: sdmt
Confirmation over: 84 (-7, +730.5) days
Baseline: fixed
Relapse influence (baseline): [30, 0] days
Relapse influence (event): [0, 0] days
Relapse influence (confirmation): [30, 0] days
Events detected: firstCDW
---
Total subjects: 7
---
Subjects with CDW: 1

Note that, in the SDMT example, visit_day is used as date column. For the MSprog() function to interpret it correctly, we need to specify date_format='day' (note that this also applies to relapse dates).

The function prints concise info (the verbose argument can be used to control the amount of info printed out – see Printing progress info section below), and generates the following two DataFrames.


1. Detailed info on each event for all subjects:

print(results_edss)
   id  nevent event_type  total_fu  time2event  bl2event  sust_days sust_last
0   1       1        CDW       534         292     292.0      242.0      True
1   2       1        CDW       730         198     198.0       84.0     False
2   3       0                  491         491       NaN        NaN       NaN
3   4       0                  586         586       NaN        NaN       NaN
4   5       1        CDW       637         140     140.0      497.0      True
5   6       0                  491         491       NaN        NaN       NaN
6   7       1        CDW       779         372     372.0      407.0      True

where: nevent is the cumulative event count for each subject; event_type characterises the event; total_fu is the number of days from start to end of follow-up; time2event is the number of days from start of follow-up to event(*); bl2event is the number of days from current baseline to event; sust_days is the number of days for which the event was sustained; sust_last reports whether the event was sustained until the last visit.

(*): For subjects with no confirmed events, time2event is the total follow-up length. To omit these subjects, set include_stable=False in MSprog().


2. Event count for each subject:

print(summary_edss)
  event_sequence  CDW
1            CDW    1
2            CDW    1
3                   0
4                   0
5            CDW    1
6                   0
7            CDW    1

Note. In this example, we called MSprog() without specifying the event argument - which defaults to firstCDW (only detect the first CDW event). For this reason,

  1. in results_edss, nevent can only be 0 or 1.

  2. in summary_edss, CDW count can only be 0 or 1, and CDI column is omitted.

In a multiple-event setting (event='multiple'), more than one event per subject can be detected, and the event-count DataFrame also includes the event sequence for each subject. See the relevant section.


Several qualitative and quantitative options for analysing disability evolution are given as optional arguments of MSprog() that can be set by the user. In order to ensure reproducibility, the results should always be complemented by the set of criteria used to obtain them. In the following sections we will go into more detail about usage and best practices for the different options. Please refer to the documentation for a complete illustration of each of the function arguments and their default values.

Outcome

In MSprog(), outcome type is specified via the mandatory argument outcome. This triggers:

  • Value range checks (e.g., EDSS in the 0-10 range)

  • Direction of worsening (e.g., for EDSS, “higher values = worse”)

  • Outcome-specific “minimum clinically relevant change”.

Outcome values are scanned in chronological order, and tested for their difference from the current reference value. Such difference is typically required to be larger than a clinically meaningful threshold \(\delta\), depending on the reference value \(x\). Default settings are implemented in pymsprog for the most used disability scales [Bosma et al., 2010, Kalinowski et al., 2022, Lorscheider et al., 2016, Strober et al., 2019]:

  • outcome='edss': Expanded Disability Status Scale (EDSS). \(\delta(x)=\begin{cases} 1.5 \quad \text{ if } x=0\\1 \quad\;\;\; \text{ if } 0 < x \leq 5\\0.5 \quad \text{ if } x>5\end{cases}\);

  • outcome='nhpt': Nine-Hole Peg Test (NHPT), for either the dominant or the non-dominant hand. \(\delta=\) 20% of reference score;

  • outcome='t25fw': Timed 25-Foot Walk (T25FW). \(\delta=\) 20% of reference score;

  • outcome='sdmt': Symbol Digit Modalities Test (SDMT). \(\delta=\) either 4 points or 20% of reference score.

These default options are internally implemented by the compute_delta() function. For example, if the baseline EDSS score is 4, a valid worsening will correspond to an increase by (at least):

from pymsprog import compute_delta
print(compute_delta(4)) # default outcome measure is 'edss'
1.0

If the baseline T25FW score is 10, the minimum shift that is considered as a valid change will be:

print(compute_delta(10, outcome='t25fw'))
2.0

Custom threshold values

The user can provide their own function to customise the computation of the threshold values, using the delta_fun argument of MSprog(). The provided function must take the baseline value as argument, and return the corresponding threshold. This applies to two scenarios.


1. The outcome of interest is among the ones listed above, but we want to define the thresholds differently. For example, we want to change the minimum \(\delta\) for SDMT to, say, “either 3 points or 10% of the reference value”. In this case, we would define:

def my_sdmt_delta(x):
    return min(3, x/10)
print('- CUSTOM minimum valid change from baseline SDMT=50: ', my_sdmt_delta(50)) # my delta
print('- DEFAULT minimum valid change from baseline SDMT=50: ', compute_delta(50, outcome='sdmt')) # default delta
- CUSTOM minimum valid change from baseline SDMT=50:  3
- DEFAULT minimum valid change from baseline SDMT=50:  4

To use our custom function, we can then set delta_fun=my_sdmt_delta in MSprog() when computing the disability events.


2. The outcome of interest is not among the ones listed above. In this case, we set outcome=None, and delta_fun as our desired custom-defined function. By default, higher values of a custom outcome are interpreted as a worsening. If the opposite direction of worsening is desired, set the worsening argument to 'decrease'. Note that the worsening argument is only used by MSprog() when outcome is set to None. Otherwise, worsening is automatically set to 'increase' if outcome is set to 'edss', 'nhpt', 't25fw', and to 'decrease' if outcome is set to 'sdmt'.

Baseline scheme

The assessment of the disability evolution strongly depends on the choices made in defining the starting point, i.e., the baseline. Different behaviours are appropriate in different contexts. This aspect is controlled by the baseline argument in MSprog(). The following alternative baseline schemes can be adopted.

  • Fixed baseline (baseline='fixed', default): the baseline visit is set to be the first available assessment.

  • Roving baseline (baseline='roving'): the baseline visit is initially set as the first available assessment, then updated after each event. This scheme is recommended in a “multiple events” setting [Kappos et al., 2018] (see example below), but not recommended for clinical trial data as it may break the randomisation.

  • Improvement-based roving baseline (baseline='roving_impr'): the baseline visit is initially set as the first available assessment, then updated after each confirmed improvement. This scheme is suitable for a first-CDW setting to discard fluctuations around the baseline [Müller et al., 2023], but not recommended for clinical trial data as it may break the randomisation.

  • Worsening-based roving baseline (baseline='roving_wors'): the baseline visit is initially set as the first available assessment, then updated after each confirmed worsening. This scheme is suitable when the endpoint of interest is a specific type of CDW: for instance, when searching specifically for PIRA events, the reference should be moved after any previous RAW event; when searching specifically for RAW events, the reference should be moved after any previous PIRA event.

NOTE. If the baseline data is stored in a separate file, the user should merge it with the main data frame containing longitudinal visit data. This can be done by inserting the baseline date and outcome value for each subject at the beginning of the data frame (not necessarily next to the other visits from the same subject – they will be grouped by subject ID).

Additional options

  • For roving baseline schemes, the re-baseline procedure can be made finer through the sub_threshold_rebl argument in MSprog(): setting sub_threshold_rebl='event', for instance, rebaseline can be triggered by any confirmed change, even if below the clinically meaningful threshold. In this case, the sub-threshold improvement events are used to move the baseline, but not listed in the event sequence.

  • When a roving baseline is used, the reference is moved after a previously detected event. The new reference may be set at the date of the previous event (proceed_from='event'), or at the first confirmation visit of such event (proceed_from='firstconf', default).

  • If the baseline visit occurs in the vicinity of a clinical relapse, the reference value may be overestimated. The relapse_to_bl argument allows to specify the minimum accepted distance (in days) of the baseline visit from the onset of a previous relapse. For instance, if baseline='fixed' and relapse_to_bl=30, the baseline visit is chosen as the first available visit satisfying the requirement of “no relapses within the preceding 30 days”. If two values are provided, they are interpreted as intervals before and after the event – e.g., relapse_to_bl=[30, 2] implements the constraint of “no relapses within the preceding 30 days or within the following 2 days”. If relapse end dates are available (may be provided as an additional column in the relapse file whose name is specified by argument renddate_col in MSprog()), the first value of relapse_to_bl is overwritten by the relapse duration, unless it was set to 0 (in which case it stays 0 and the baseline is not moved based on the influence of previous relapses). If the baseline argument is set to a roving scheme, all the above applies to every baseline change as well.

  • On top of the chosen baseline scheme, post-relapse re-baseline [Kappos et al., 2020] can be applied by setting relapse_rebl=True in MSprog(). If this is enabled, the onset of each relapse prompts a re-baseline to the next available visit following the relapse and out of its influence (as per relapse_to_bl). As an additional constraint, an updated baseline can be forced to have a score greater than or equal to the previous baseline by setting bl_geq=True in MSprog().

  • Due to fluctuations in the data, the new reference value following a re-baseline may fall at a local minimum or maximum of the disability trajectory. These may be excluded for consideration as baseline visits using the skip_local_extrema argument in MSprog().


As already mentioned above, extra caution should be used when applying any re-baseline scheme to randomised data, as moving the reference value based on post-randomisation events can introduce bias (especially if the occurrence of these events is influenced by the treatment). For clinical trial data, general recommendation is to use a fixed baseline (at the time of randomisation).

Multiple events

The event argument allows to specify which events to detect. By default, it is set to 'firstCDW' (only detect the first CDW event). It can be set to 'multiple' to sequentially detect all CDW or CDI events.

For example, extracting multiple EDSS events for subject 4 from toydata_visits with a fixed baseline would result in the following.

print('\nData:')
print(toydata_visits.loc[toydata_visits['id']==4, ['date', 'EDSS']].reset_index(drop=True)) # EDSS visits

_, results = MSprog(toydata_visits, 'id', 'EDSS', 'date', outcome='edss',
                    relapse=toydata_relapses, 
                    subjects=[4],
                    conf_days=12 * 7, 
                    event='multiple', baseline='fixed',  # <---
                    verbose=0)
print('\n---Results with fixed baseline:---')
print(results.drop(columns='id').set_index('nevent', drop=True).T) # results
Data:
         date  EDSS
0  2021-09-18   4.5
1  2021-12-04   3.5
2  2022-03-12   3.5
3  2022-07-19   5.0
4  2022-10-05   5.0
5  2023-01-16   5.5
6  2023-04-27   5.0

---Results with fixed baseline:---
nevent          1
event_type    CDI
total_fu      586
time2event     77
bl2event     77.0
sust_days      98
sust_last   False

With these settings, the EDSS improvement at visit 2 (EDSS=3.5, confirmed at visit 3) does not trigger a re-baseline. The algorithm keeps searching for events from after the confirmation visit, evaluating the changes relative to the original baseline (EDSS=4.5). The subsequent EDSS worsening is therefore not detected. On the other hand, adopting a roving baseline scheme we get:

_, results = MSprog(toydata_visits, 'id', 'EDSS', 'date', outcome='edss',
                    relapse=toydata_relapses, 
                    subjects=[4],
                    conf_days=12 * 7, 
                    event='multiple', baseline='roving',  # <-----
                    verbose=0)
print('\n---Results with roving baseline:---')
print(results.drop(columns='id').set_index('nevent', drop=True).T) # results
---Results with roving baseline:---
nevent          1      2
event_type    CDI    CDW
total_fu      586    586
time2event     77    304
bl2event     77.0  129.0
sust_days      98    282
sust_last   False   True

The baseline has now been moved to the confirmation visit of the first CDI event (visit 3, EDSS=3.5). With respect to this baseline, a subsequent CDW (PIRA) event is detected at visit 4 (EDSS=5), confirmed at visit 5.

Other valid values for the event argument are:

  • 'first': only the very first confirmed event – CDI or CDW

  • 'firstCDI': first CDI event

  • 'firstPIRA': first PIRA event

  • 'firstRAW': first RAW event.

Note: if, for example, 'event=firstPIRA', all non-PIRA events preceding the first PIRA are actually detected – although not reported. This allows to combine this scenario with a roving baseline scheme, where the baseline is moved after every confirmed non-PIRA event. For instance, if a RAW event occurs first, the subsequent worsening would be counted from the re-baselined outcome value after the RAW.

Event confirmation

An event is only validated if the change in the outcome value from the current baseline is maintained up to(*) a subsequent confirmation visit at a pre-specified distance from the event [Ontaneda et al., 2017]. The event is confirmed if the difference in the outcome value from the baseline score remains above-threshold at all assessments up to the confirmation visit. For example, with the default EDSS thresholds, an increase in EDSS from 4 points to 6 points is confirmed if EDSS=5 at all subsequent visits; it is not confirmed if EDSS=4.5 at all subsequent visits.

The chosen confirmation period depends on the type of study and on the frequency of visits, and can be set in MSprog() by using the conf_days argument. A tolerance interval around conf_days can be specified using the conf_tol_days argument, given as a sequence of two integers (left and right tolerance)(**). Setting the right end of the interval to Inf allows event confirmation to occur at any visit after conf_days. Let’s see two examples.

i. A common setting for clinical trial data would be: conf_days=7*12, conf_tol_days=[0, np.inf], i.e., “confirmation over 12 or more weeks”.

ii. A common setting for observational data would be: conf_days=7*24, conf_tol_days=[0, 7*12], i.e., the confirmation visit must lie between 24 weeks after the event, and 36 weeks after the event.

Note. “Confirmed over 12 or more weeks” (i.e., no upper bound) \(\neq\) “sustained over the whole follow up”.
Setting no upper bound to the confirmation interval means that the first visit at least 12 weeks after the initial change is selected as confirmation visit. The confirmation window is then defined as the interval between the initial worsening and that confirmation visit, and scores are checked therein.

Additional options

  • If conf_days is specified as a sequence of values (e.g., conf_days=[7*12, 7*24]), events are retained if confirmed by at least one visit falling within any of the specified periods (here, 12 or 24 weeks with their relative tolerance interval)(***). Confirmation dates/values at all specified values will be included in the results table (if include_dates=True / include_values=True in MSprog()).

  • If the confirmation visit occurs in the vicinity of a clinical relapse, it is typically considered invalid. The relapse_to_conf argument allows to specify the minimum accepted distance (in days) of a confirmation visit from the onset of a previous relapse. For example, relapse_to_conf=30 implements the constraint “no relapses within the preceding 30 days” for a visit to be used as confirmation. If two values are provided, they are interpreted as intervals before and after the event – e.g., relapse_to_conf=[30, 2] implements the constraint “no relapses within the preceding 30 days or within the following 2 days”. If relapse end dates are available (may be provided as an additional column in the relapse file whose name is specified by argument renddate_col in MSprog()), the first value of relapse_to_conf is overwritten by the relapse duration, unless it was set to 0 (in which case it stays 0 and confirmation visits are not discarded based on the influence of previous relapses).

  • By default, a disability worsening occurring at the last available assessment is ignored by MSprog(). This behaviour can be changed using the impute_last_visit argument. If set to TRUE, any disability worsening occurring at the last visit is retained even though unconfirmed. If set to a numeric value N, unconfirmed worsening events are included only if occurring within N days of follow up (e.g., in case of early discontinuation).

  • In scenarios (e.g., a clinical trial) where disability course is compared between two cohorts with different relapse rates, there may be an imbalance in visit frequency between the two groups, due to unscheduled visits after relapses. This implies a higher probability of detecting outcome changes in the group with more assessments. Such issue is sometimes addressed by only using scheduled visits as confirmation visits [Hauser et al., 2017]. It is possible to implement this with MSprog() by including an additional column in the longitudinal data, specifying which visits can (True) or cannot (False) be used as confirmation visits. The name of such column must be specified using the validconf_col argument – for example, if the database contains a column 'scheduled' tracking scheduled visits, unscheduled visits can be excluded from consideration as confirmation visits by setting validconf_col='scheduled'in MSprog().

(*) The value change from baseline must also be maintained at all intermediate visits between the event and the confirmation visit. This behaviour may be changed by setting the check_intermediate argument to False (in this case, the change only needs to be confirmed at the designated confirmation visit). We do not recommend this, as it may lead to including random fluctuations as events. We included this option to provide the possibility of replicating the results from previous studies that used this rationale.

(**) conf_tol_days can also be specified as a single integer if the same tolerance on left and right is desired.

(***) An event is only confirmed if the value change from baseline is maintained at all visits up to the confirmation visit. So an event can only be confirmed over 24 weeks but not over 12 weeks if there are no valid confirmation visits falling within the 12-week window (unless check_intermediate=FALSE).

Sustained CDW or CDI

In addition to event confirmation, some studies require events to be sustained for a specified period of time. The require_sust_days argument allows to specify, if desired, the length of such period. For example, if require_sust_days=7*48, an event is only retained if the change in the outcome measure from the current baseline is confirmed at all visits falling within the following 48 weeks. If the event is sustained for the entire follow-up, it is retained even if the follow-up period ends <48 weeks after the event. Setting require_sust_days=Inf, events are retained only when sustained for the entire follow-up duration.

The require_sust_days argument may be of use in the following scenarios.

  • Suppose the maximum follow-up in the population is 96 weeks excluding the baseline visit. If one wants to detect CDWs sustained over the whole follow-up, coding it as “worsening confirmed over 96 weeks or more” (conf_days=7*96, conf_tol_days=c(0, Inf)) would discard CDWs of patients whose follow-up is shorter than 96 weeks. In a scenario where patients have different follow-up lengths, the require_sust_days argument offers a more reliable way of implementing this, e.g., as “worsening confirmed over 12 weeks or more and sustained for the whole follow-up period” (conf_days=7*12, conf_tol_days=c(0, Inf), require_sust_days=Inf).

  • In observational data with coarsely- and unevenly-spaced visits, the require_sust_days argument may be used to detect sustained CDW (e.g., require_sust_days=7*96) while still requiring the presence of a confirmation visit within a specified interval from the event (e.g., conf_days=7*12, conf_tol_days=c(0,7*12)).

  • Sustained disability worsening can be required to exclude transient RAW events.

Relapse-based classification of CDW events

Detected CDW events may be further categorised based on their timing with respect to relapses (*) and identified as RAW or PIRA. This has to be explicitly enabled by setting the flag RAW_PIRA=True in MSprog() (unless the endpoint of interest is already firstPIRA or firstRAW). On top of that, the definitions used to label a CDW as RAW or as PIRA are controlled by arguments relapse_assoc and relapse_indep, as detailed below.

(*) Relapse data is to be provided using the optional relapse argument in MSprog(). It should be given as a data frame containing subject IDs and relapse onset dates (see section Input data). If the names of columns with subject IDs and dates in the relapse database are different from the main database, they must be specified using arguments rsubj_col and rdate_col.

RAW

RAW events are typically defined as CDW events occurring within a specified interval from the onset of a relapse. The length (in days) of such interval can be specified using the relapse_assoc argument in MSprog(). Common values are 30 or 90 days. Additionally, one may also provide a maximum distance from relapses whose onset is after the event. The logic is that the reported onset date of a relapse may be slightly delayed. The examples below illustrate the usage of the relapse_assoc argument:

  • relapse_assoc=90 (equivalent to relapse_assoc=[90, 0]): RAW = “a CDW event occurring within 90 days after the onset of a relapse”;

  • relapse_assoc=[90, 7]: RAW = “a CDW event occurring no more than 90 days after and no more than 7 days before the onset of a relapse”.

PIRA

In the literature, PIRA is typically defined by requiring an absence of relapses within appropriate intervals anchored to the dates of the baseline visit, the worsening event, and the confirmation visit [Cagol et al., 2022, Kappos et al., 2020, Müller et al., 2023].

The relapse_indep argument in MSprog() allows to specify custom relapse-free intervals. The argument must be provided in the form of a dictionary as follows:

{'prec': (p0, p1), 'event': (e0, e1), 'conf': (c0, c1)}.

The dictionary specifies the intervals [p0, p1], [e0, e1], and [c0, c1] around (any subset of) the three checkpoints:

  • 'prec': a visit preceding the event – see below;

  • 'event': the disability worsening event;

  • 'conf': the first available confirmation visit.

The dictionary can also optionally contain a key-value pair specifying how to choose 'prec':

  • 'prec_type': 'baseline''prec' is the current baseline;

  • 'prec_type': 'last''prec' is the last visit before the event;

  • 'prec_type': 'last_lower''prec' is the last visit before event onset with a clinically meaningful score difference from it: i such that outcome[event] - outcome[i] >= delta_fun(outcome[i]) (if worsening='increase', the opposite otherwise) and same for the confirmation visit.

If 'prec_type' is not in the dictionary keys, it will be automatically assigned as 'baseline'.

If both ends of an interval are 0 (e.g., if both p0=0 and p1=0), the checkpoint is ignored. To merge two intervals together, set both the right end of the first interval and the left end of the second interval to None (e.g., “between baseline and event onset”: p1=None and e0=None).

For example, in Müller et al. [2023], the authors recommend an absence of relapses during the 90 days before and 30 days after the event, and during the 90 days before and 30 days after the confirmation visit for a CDW event to be considered as PIRA. This translates into: p0 = 0, p1 = 0, e0 = 90, e1 = 30, c0 = 90, c1 = 30. When a high specificity is desired, they recommend an absence of relapses in the whole period between the reference visit and the confirmation visit, which is implemented by setting: p0 = 0, p1 = None, e0 = None, e1 = None, c0 = None, c1 = 0.

Example

The following code detects multiple events with the RAW_PIRA flag enabled.

summary, results = MSprog(data=toydata_visits,
                 subj_col='id', value_col='EDSS', date_col='date', outcome='edss', 
                 event='multiple', RAW_PIRA=True, baseline='roving', 
                 relapse=toydata_relapses) 
---
Outcome: edss
Confirmation over: 84 (-7, +730.5) days
Baseline: roving
Relapse influence (baseline): [30, 0] days
Relapse influence (event): [0, 0] days
Relapse influence (confirmation): [30, 0] days
Events detected: multiple
---
Total subjects: 7
---
Subjects with CDW: 5 (PIRA: 5; RAW: 1)
Subjects with CDI: 2
---
CDW events: 6 (PIRA: 5; RAW: 1)
CDI events: 2

The summary DataFrame (and the info printed out by the function) now includes RAW and PIRA counts.

print(summary)
  event_sequence  CDI  CDW  RAW  PIRA
1           PIRA    0    1    0     1
2      RAW, PIRA    0    2    1     1
3                   0    0    0     0
4      CDI, PIRA    1    1    0     1
5           PIRA    0    1    0     1
6            CDI    1    0    0     0
7           PIRA    0    1    0     1

In the results DataFrame, each CDW event is now further characterised by a CDW_type column:

print(results)
   id  nevent event_type CDW_type  total_fu  time2event  bl2event  sust_days  \
0   1       1        CDW     PIRA       534         292     292.0      242.0   
1   2       1        CDW      RAW       730         198     198.0       84.0   
2   2       2        CDW     PIRA       730         539     257.0      191.0   
3   3       0                           491         491       NaN        NaN   
4   4       1        CDI                586          77      77.0       98.0   
5   4       2        CDW     PIRA       586         304     129.0      282.0   
6   5       1        CDW     PIRA       637         140     140.0      497.0   
7   6       1        CDI                491         120     120.0      232.0   
8   7       1        CDW     PIRA       779         372     372.0      407.0   

  sust_last  
0      True  
1     False  
2      True  
3       NaN  
4     False  
5      True  
6      True  
7     False  
8      True  

MSprog() outputs

What to include in results

The results DataFrame provides extended info on each event for all subjects. The following arguments can be used to control the information included.

  • include_dates: if True, the results will include the dates of: event onset; baseline; last visit before event onset at clinically meaningful score difference from it; confirmation visit.

  • include_value: if True, the results will include the outcome value at: event onset; baseline; last visit before event onset at clinically meaningful score difference from it; confirmation visit.

  • include_stable: if True, subjects with no confirmed events are included in the results, with time2event = total follow up.

For example:

summary, results = MSprog(toydata_visits,
                         subj_col='id', value_col='EDSS', date_col='date',
                         outcome='edss',
                         relapse=toydata_relapses,
                         include_dates=True, include_value=True, include_stable=False,
                         verbose=0)
print(results.T)
                                    0                    1  \
id                                  1                    2   
nevent                              1                    1   
event_type                        CDW                  CDW   
date              2022-07-12 00:00:00  2021-06-12 00:00:00   
value                             6.0                  5.5   
bl_date           2021-09-23 00:00:00  2020-11-26 00:00:00   
bl_value                          5.5                  4.0   
last_delta_date   2022-04-27 00:00:00  2020-12-30 00:00:00   
last_delta_value                  5.5                  4.0   
total_fu                          534                  730   
time2event                        292                  198   
bl2event                        292.0                198.0   
conf84_date       2022-11-06 00:00:00  2021-09-04 00:00:00   
conf84_value                      6.0                  5.0   
sust_days                         242                   84   
sust_last                        True                False   

                                    2                    3  
id                                  5                    7  
nevent                              1                    1  
event_type                        CDW                  CDW  
date              2021-11-01 00:00:00  2022-03-17 00:00:00  
value                             5.0                  4.5  
bl_date           2021-06-14 00:00:00  2021-03-10 00:00:00  
bl_value                          4.0                  3.5  
last_delta_date   2021-09-12 00:00:00  2021-12-02 00:00:00  
last_delta_value                  4.0                  3.5  
total_fu                          637                  779  
time2event                        140                  372  
bl2event                        140.0                372.0  
conf84_date       2022-01-17 00:00:00  2022-07-22 00:00:00  
conf84_value                      5.5                  4.5  
sust_days                         497                  407  
sust_last                        True                 True  

Printing progress info

The verbose argument controls the amount of info printed out by the MSprog() function. When setting verbose=0, the function prints no info. When setting verbose=1 (default), the function only prints out concise info. When setting verbose=2, the function prints out an extended log of the ongoing computations.

For further insight into the event detection process, the return_unconfirmed argument allows to also return all score changes from baseline that were not identified as valid events (e.g., not confirmed).

See the example below.

summary, results, unconfirmed = MSprog(toydata_visits,
                          subj_col='id', value_col='EDSS', date_col='date',
                          outcome='edss',
                          event='multiple', baseline='roving',
                          relapse=toydata_relapses,
                          return_unconfirmed=True,  # <---
                          verbose=2)  # <---
Subject #1: 8 visits, 0 relapses
Baseline at visit no.1
Searching for events from visit no.2 on
EDSS change at visit no.5 (2022-07-12); potential confirmation visits available: no.6, 7, 8
Found EDSS-CDW (visit no.5, 2022-07-12) confirmed at 84 weeks, sustained up to visit no.8 (2023-03-11)
Baseline at visit no.6
Searching for events from visit no.7 on
No EDSS change in any subsequent visit: end process
Event sequence: CDW

Subject #2: 10 visits, 2 relapses
Baseline at visit no.1
Searching for events from visit no.2 on
EDSS change at visit no.4 (2021-06-12); potential confirmation visits available: no.5, 6, 7, 8, 9, 10
Found EDSS-CDW (visit no.4, 2021-06-12) confirmed at 84 weeks, sustained up to visit no.5 (2021-09-04)
Baseline at visit no.5
Searching for events from visit no.6 on
EDSS change at visit no.8 (2022-05-19); potential confirmation visits available: no.9, 10
Found EDSS-CDW (visit no.8, 2022-05-19) confirmed at 84 weeks, sustained up to visit no.10 (2022-11-26)
Baseline at visit no.9
Searching for events from visit no.10 on
No EDSS change in any subsequent visit: end process
Event sequence: CDW, CDW
Subject #3: 7 visits, 1 relapse
Baseline at visit no.1
Searching for events from visit no.2 on
EDSS change at visit no.2 (2021-12-01); potential confirmation visits available: no.4, 5, 7
Change not confirmed: proceed with search
Searching for events from visit no.3 on
No EDSS change in any subsequent visit: end process
Event sequence: -
Subject #4: 7 visits, 0 relapses
Baseline at visit no.1
Searching for events from visit no.2 on
EDSS change at visit no.2 (2021-12-04); potential confirmation visits available: no.3, 4, 5, 6, 7
Found EDSS-CDI (visit no.2, 2021-12-04) confirmed at 84 weeks, sustained up to visit no.3 (2022-03-12)
Baseline at visit no.3
Searching for events from visit no.4 on
EDSS change at visit no.4 (2022-07-19); potential confirmation visits available: no.5, 6, 7
Found EDSS-CDW (visit no.4, 2022-07-19) confirmed at 84 weeks, sustained up to visit no.7 (2023-04-27)
Baseline at visit no.5
Searching for events from visit no.6 on
No EDSS change in any subsequent visit: end process
Event sequence: CDI, CDW

Subject #5: 9 visits, 0 relapses
Baseline at visit no.1
Searching for events from visit no.2 on
EDSS change at visit no.3 (2021-11-01); potential confirmation visits available: no.4, 5, 6, 7, 8, 9
Found EDSS-CDW (visit no.3, 2021-11-01) confirmed at 84 weeks, sustained up to visit no.9 (2023-03-13)
Baseline at visit no.4
Searching for events from visit no.5 on
EDSS change at visit no.5 (2022-04-30); potential confirmation visits available: no.7, 8, 9
Change not confirmed: proceed with search
Searching for events from visit no.6 on
EDSS change at visit no.8 (2022-12-15); potential confirmation visits available: no.9
Change not confirmed: proceed with search
Searching for events from visit no.9 on
No EDSS change in any subsequent visit: end process
Event sequence: CDW

Subject #6: 7 visits, 1 relapse
Baseline at visit no.1
Searching for events from visit no.2 on
EDSS change at visit no.3 (2022-02-15); potential confirmation visits available: no.4, 5, 7
Found EDSS-CDI (visit no.3, 2022-02-15) confirmed at 84 weeks, sustained up to visit no.5 (2022-10-05)
Baseline at visit no.4
Searching for events from visit no.5 on
EDSS change at visit no.6 (2022-12-22); potential confirmation visits available: no.
Change not confirmed: proceed with search
Searching for events from visit no.7 on
EDSS change at visit no.7 (2023-02-21); potential confirmation visits available: no.
Change not confirmed: proceed with search
Searching for events from visit no.- on
No EDSS change in any subsequent visit: end process
Event sequence: CDI

Subject #7: 9 visits, 1 relapse
Baseline at visit no.1
Searching for events from visit no.2 on
EDSS change at visit no.3 (2021-09-11); potential confirmation visits available: no.4, 5, 6, 7, 8, 9
Change not confirmed: proceed with search
Searching for events from visit no.4 on
EDSS change at visit no.5 (2022-03-17); potential confirmation visits available: no.6, 7, 8, 9
Found EDSS-CDW (visit no.5, 2022-03-17) confirmed at 84 weeks, sustained up to visit no.9 (2023-04-28)
Baseline at visit no.6
Searching for events from visit no.7 on
No EDSS change in any subsequent visit: end process
Event sequence: CDW

---
Outcome: edss
Confirmation over: 84 (-7, +730.5) days
Baseline: roving
Relapse influence (baseline): [30, 0] days
Relapse influence (event): [0, 0] days
Relapse influence (confirmation): [30, 0] days
Events detected: multiple
---
Total subjects: 7
---
Subjects with CDW: 5
Subjects with CDI: 2
---
CDW events: 6
CDI events: 2

Let’s visualise the unconfirmed events:

print('---Unconfirmed events:---')
print(unconfirmed)
---Unconfirmed events:---
   id       date  value    bl_date  bl_value  closest_rel-  closest_rel+
0   3 2021-12-01    4.5 2021-10-18       3.5           inf         365.0
1   5 2022-04-30    5.0 2022-01-17       5.5           inf           inf
2   5 2022-12-15    5.0 2022-01-17       5.5           inf           inf
3   6 2022-12-22    6.0 2022-09-03       4.5           4.0           inf
4   6 2023-02-21    6.0 2022-09-03       4.5          65.0           inf
5   7 2021-09-11    5.0 2021-03-10       3.5          -0.0           0.0

Time to disability milestone

Instead of studying disability course with respect to a baseline value, one can focus on the time taken to reach a specific disability milestone (e.g., EDSS \(\geq\) 6.0). This can be computed using function pymsprog.value_milestone().

The following code detects the time to EDSS \(\geq\) 6.0 for each subject in the toy data.

from pymsprog import value_milestone

vm = value_milestone(toydata_visits, milestone=6,
                 subj_col='id', value_col='EDSS', date_col='date', outcome='edss',
                 relapse=toydata_relapses,
                 verbose=0)

The function returns a DataFrame indexed by subject IDs:

print(vm)
        date  EDSS  time2event  observed
1 2022-07-12   6.0       292.0      True
2 2022-05-19   6.0       539.0      True
3 2023-02-21   NaN       491.0     False
4 2023-04-27   NaN       586.0     False
5 2023-03-13   NaN       637.0     False
6 2023-02-21   NaN       491.0     False
7 2023-04-28   NaN       779.0     False

where: 'date' contains the date of the first confirmed EDSS \(\geq\) 6.0 (or last date of follow-up if milestone is not reached or not confirmed); 'EDSS' contains the first EDSS value \(\geq\) 6.0, if present, otherwise no value; 'time2event' contains the time to reach EDSS \(\geq\) 6.0 (or total follow-up length if not reached or not confirmed); 'observed' indicates whether EDSS \(\geq\) 6.0 was reached (1) or not (0).

Several arguments controlling the behaviour of the MSprog() function are also available for the value_milestone() function (e.g., require_sust_days, impute_last_visit…). Please refer to the documentation for a complete illustration of each of the function arguments and their default values.

References

[1]

L V A E Bosma, J J Kragt, L Brieva, Z Khaleeli, X Montalban, C H Polman, A J Thompson, M Tintoré, and B M J Uitdehaag. Progression on the multiple sclerosis functional composite in multiple sclerosis: what is the optimal cut-off for the three components? Mult Scler, 16(7):862–867, Jul 2010. doi:10.1177/1352458510370464.

[2]

Alessandro Cagol, Sabine Schaedelin, Muhamed Barakovic, Pascal Benkert, Ramona-Alexandra Todea, Reza Rahmanzadeh, Riccardo Galbusera, Po-Jui Lu, Matthias Weigel, Lester Melie-Garcia, Esther Ruberte, Nina Siebenborn, Marco Battaglini, Ernst-Wilhelm Radue, Özgür Yaldizli, Johanna Oechtering, Tim Sinnecker, Johannes Lorscheider, Bettina Fischer-Barnicol, Stefanie Müller, Lutz Achtnichts, Jochen Vehoff, Giulio Disanto, Oliver Findling, Andrew Chan, Anke Salmen, Caroline Pot, Claire Bridel, Chiara Zecca, Tobias Derfuss, Johanna M. Lieb, Luca Remonda, Franca Wagner, Maria I. Vargas, Renaud Du Pasquier, Patrice H. Lalive, Emanuele Pravatà, Johannes Weber, Philippe C. Cattin, Claudio Gobbi, David Leppert, Ludwig Kappos, Jens Kuhle, and Cristina Granziera. Association of Brain Atrophy With Disease Progression Independent of Relapse Activity in Patients With Relapsing Multiple Sclerosis. JAMA Neurology, 79(7):682–692, 07 2022. URL: https://doi.org/10.1001/jamaneurol.2022.1025, arXiv:https://jamanetwork.com/journals/jamaneurology/articlepdf/2792415/jamaneurology\_cagol\_2022\_oi\_220020\_1656444982.99513.pdf, doi:10.1001/jamaneurol.2022.1025.

[3]

Stephen L Hauser, Amit Bar-Or, Giancarlo Comi, Gavin Giovannoni, Hans-Peter Hartung, Bernhard Hemmer, Fred Lublin, Xavier Montalban, Kottil W Rammohan, Krzysztof Selmaj, Anthony Traboulsee, Jerry S Wolinsky, Douglas L Arnold, Gaelle Klingelschmitt, Donna Masterman, Paulo Fontoura, Shibeshih Belachew, Peter Chin, Nicole Mairon, Hideki Garren, and Ludwig Kappos. Ocrelizumab versus interferon beta-1a in relapsing multiple sclerosis. N Engl J Med, 376(3):221–234, Jan 2017. doi:10.1056/NEJMoa1601277.

[4]

Anissa Kalinowski, Gary Cutter, Nina Bozinov, Jessica A Hinman, Michael Hittle, Robert Motl, Michelle Odden, and Lorene M Nelson. The timed 25-foot walk in a large cohort of multiple sclerosis patients. Mult Scler, 28(2):289–299, Feb 2022. doi:10.1177/13524585211017013.

[5]

Ludwig Kappos, Helmut Butzkueven, Heinz Wiendl, Timothy Spelman, Fabio Pellegrini, Yi Chen, Qunming Dong, Harold Koendgen, Shibeshih Belachew, Maria Trojano, and On Behalf of the Tysabri® Observational Program (TOP) Investigators. Greater sensitivity to multiple sclerosis disability worsening and progression events using a roving versus a fixed reference value in a prospective cohort study. Multiple Sclerosis Journal, 24(7):963–973, 2018. PMID: 28554238. URL: https://doi.org/10.1177/1352458517709619, arXiv:https://doi.org/10.1177/1352458517709619, doi:10.1177/1352458517709619.

[6] (1,2,3)

Ludwig Kappos, Jerry S Wolinsky, Gavin Giovannoni, Douglas L Arnold, Qing Wang, Corrado Bernasconi, Fabian Model, Harold Koendgen, Marianna Manfrini, Shibeshih Belachew, and Stephen L Hauser. Contribution of relapse-independent progression vs relapse-associated worsening to overall confirmed disability accumulation in typical relapsing multiple sclerosis in a pooled analysis of 2 randomized clinical trials. JAMA Neurol, 77(9):1132–1140, Sep 2020. doi:10.1001/jamaneurol.2020.1568.

[7]

Johannes Lorscheider, Katherine Buzzard, Vilija Jokubaitis, Tim Spelman, Eva Havrdova, Dana Horakova, Maria Trojano, Guillermo Izquierdo, Marc Girard, Pierre Duquette, Alexandre Prat, Alessandra Lugaresi, François Grand'Maison, Pierre Grammond, Raymond Hupperts, Raed Alroughani, Patrizia Sola, Cavit Boz, Eugenio Pucci, Jeanette Lechner-Scott, Roberto Bergamaschi, Celia Oreja-Guevara, Gerardo Iuliano, Vincent Van Pesch, Franco Granella, Cristina Ramo-Tello, Daniele Spitaleri, Thor Petersen, Mark Slee, Freek Verheul, Radek Ampapa, Maria Pia Amato, Pamela McCombe, Steve Vucic, JoséLuis Sánchez Menoyo, Edgardo Cristiano, Michael H Barnett, Suzanne Hodgkinson, Javier Olascoaga, Maria Laura Saladino, Orla Gray, Cameron Shaw, Fraser Moore, Helmut Butzkueven, and Tomas Kalincik. Defining secondary progressive multiple sclerosis. Brain, 139(Pt 9):2395–2405, Sep 2016. doi:10.1093/brain/aww173.

[8]

Fred D. Lublin, Stephen C. Reingold, Jeffrey A. Cohen, Gary R. Cutter, Per Soelberg Sørensen, Alan J. Thompson, Jerry S. Wolinsky, Laura J. Balcer, Brenda Banwell, Frederik Barkhof, Jr Bruce Bebo, Peter A. Calabresi, Michel Clanet, Giancarlo Comi, Robert J. Fox, Mark S. Freedman, Andrew D. Goodman, Matilde Inglese, Ludwig Kappos, Bernd C. Kieseier, John A. Lincoln, Catherine Lubetzki, Aaron E. Miller, Xavier Montalban, Paul W. O\textquoteright Connor, John Petkau, Carlo Pozzilli, Richard A. Rudick, Maria Pia Sormani, Olaf Stüve, Emmanuelle Waubant, and Chris H. Polman. Defining the clinical course of multiple sclerosis. Neurology, 83(3):278–286, 2014. URL: https://n.neurology.org/content/83/3/278, arXiv:https://n.neurology.org/content/83/3/278.full.pdf, doi:10.1212/WNL.0000000000000560.

[9] (1,2,3)

Jannis Müller, Alessandro Cagol, Johannes Lorscheider, Charidimos Tsagkas, Pascal Benkert, Özgür Yaldizli, Jens Kuhle, Tobias Derfuss, Maria Pia Sormani, Alan Thompson, Cristina Granziera, and Ludwig Kappos. Harmonizing definitions for progression independent of relapse activity in multiple sclerosis: a systematic review. JAMA Neurol, 80(11):1232–1245, Nov 2023. doi:10.1001/jamaneurol.2023.3331.

[10]

Daniel Ontaneda, Alan J Thompson, Robert J Fox, and Jeffrey A Cohen. Progressive multiple sclerosis: prospects for disease therapy, repair, and restoration of function. Lancet, 389(10076):1357–1366, Apr 2017. doi:10.1016/S0140-6736(16)31320-4.

[11]

Lauren Strober, John DeLuca, Ralph Hb Benedict, Adam Jacobs, Jeffrey A Cohen, Nancy Chiaravalloti, Lynn D Hudson, Richard A Rudick, and Nicholas G LaRocca. Symbol digit modalities test: a valid clinical trial endpoint for measuring cognition in multiple sclerosis. Mult Scler, 25(13):1781–1790, Nov 2019. doi:10.1177/1352458518808204.