{'fraction_non_alphanumeric': 0.034618140134689786, 'fraction_numerical': 0.013554950565983666, 'mean_word_length': 4.589620374819798, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Feature Importance Explanations for Temporal Black-Box Models

Models in the supervised learning framework may capture rich and complex representations over the features that are hard for humans to interpret. Existing methods to explain such models are often specific to architectures and data where the features do not have a time-varying component. In this work, we propose TIME, a method to explain models that are inherently temporal in nature. Our approach (i) uses a model-agnostic permutation-based approach to analyze global feature importance, (ii) identifies the importance of salient features with respect to their temporal ordering as well as localized windows of influence, and (iii) uses hypothesis testing to provide statistical rigor.


Introduction
The last decade has seen an explosion in models that learn rich representations over large, complex parameter spaces. These have increasingly been applied in domains with a high degree of social impact, such as healthcare, but this very complexity makes them black-boxes whose decisionmaking is hard to explain, a critical deficit in many such domains. There has thus been a concomitant rise in methods to generate explanations for black-box models. Existing research has largely focused on explaining models trained over tabular data, where each feature takes a single value per instance, instead of explaining temporal models, where the instances consist of sequences or time series.
Most explanation methods, including model-agnostic feature importance methods such as LIME (Ribeiro et al., 2016) and SHAP (Lundberg & Lee, 2017), are designed for tabular representations. Ismail et al. (2020) demonstrate the unreliability and inaccuracy of several explanation methods designed for tabular data in identifying feature importance in temporal models. Some approaches have focused on interpreting specific temporal model architectures, such as recurrent neural networks (Karpathy et al., 2015;Suresh et al., 2017;Ismail et al., 2019) and attention-based models (Choi et al., 2016;Zhang et al., 2019), while others have explored methods to encourage temporal models during training to be more interpretable using tree regularization (Wu et al., 2017) and game-theoretic characterizations (Lee et al., 2018). However, model-agnostic explanation for temporal models has begun to be addressed only recently. Tonekaboni et al. (2020) propose FIT, a method to assign importance scores to observations over time for sequence-sequence models, and Bento et al. (2020) propose TimeSHAP, an extension of SHAP to temporal models. Importantly, these methods focus on local interpretability, which seeks to explain individual predictions in terms of their important features, rather than global interpretability, which seeks to identify features important to the model as a whole (Ibrahim et al., 2019;Covert et al., 2020b).
Our approach is most similar to permutation-based methods that identify global feature importance for learned models. Breiman (2001) uses permutations to identify important features in random forests, and many variants of feature importance using permutation tests have since been studied (Strobl et al., 2008;Ojala & Garriga, 2010;Altmann et al., 2010;Gregorutti et al., 2015;Fisher et al., 2019;Zhou & Hooker, 2020). The simplicity and generality of permutation tests makes them attractive as a tool for model-agnostic explanation. However, existing methods have focused on permutations of individual features, or in some cases groups of features, as part of a tabular representation. We extend permutation tests to temporal models by using the temporal structure of the data to guide the selection of features to be permuted and the resulting explanations.
In this work, we propose Temporal Importance Model Explanation (TIME), a method for explaining temporal blackbox models. Our approach is model-agnostic, produces global explanations, and elicits specific properties of temporal models. It takes as input a learned model over features representing sequences or time-series, and a test data set used to evaluate the model, and does the following: (i) it identifies features that are important for the model's predictions across the distribution of instances, (ii) for each such feature, it identifies the most important temporal window that the model focuses on, (iii) it determines whether the model's predictions are dependent on the ordering of the values within the window, (iv) it uses hypothesis testing and (a) Figure 1. (a) Time series for positive (green) and negative (red) instances for four different features, illustrating temporal properties of the features that a learned model may capture. (b) A trained binary classification model over the four time-varying features, whose underlying function uses the features' temporal properties to capture the target concept. xA is not used by the model; all timesteps for xB are equally important; the model focuses on windows [c1, c2] and [d1, d2] for xC and xD respectively; the ordering of values is important only for xD. (c) The output of TIME, showing for each feature (i) its overall importance to the model, (ii) the most important window that the model focuses on, and (iii) whether the ordering of the values within the window is important to the model. a false discovery rate control methodology to identify important features and their temporal properties with statistical rigor, and (v) it treats the model as a black-box, and thus and may be used to analyze a wide variety of temporal or sequential models. Figure 1 illustrates the setting and our approach.

Identifying Important Features/Timesteps via Permutation
Non-temporal models. We first outline the case of a model trained on a tabular data set where each feature takes a single value per instance. Consider a model f over D features, trained to predict a target y. We are interested in examining the importance of a given feature j for the model in predicting y. To focus on the model's generalization performance, we assume a test set comprising M instances is available to analyze the model. Let (x (i) , y (i) ) be the i th instance-target pair, and L be a loss function linking the model output f (x) to the target y.
We define the baseline output of the model f for instance i as: The perturbed output of the model for instance i w.r.t feature j and another instance l = i is defined as: where the value of feature j is replaced by its corresponding value from instance l, as shown in Figure 2a. Then, we can compute the change in loss between the perturbed and baseline losses as: Let S i ⊂ {1, 2, . . . M }\{i} be a subset of instances excluding i. We replace the value of feature j for instance i with its value from each instance in S i and compute the average change in loss as: By averaging this calculation over all instances i = 1 . . . M using subsets S i selected by permuting the data set, we compute the importance score of feature j as: A model includes many features, all of which may have some effect on the model's output, but only some of which are useful in predicting the target. We consider a feature to be important if it has a positive association with the target, in the sense that permuting the value of the feature increases the model's loss on average. We use hypothesis testing to capture this notion of importance, as outlined in Subsection 2.4.
Temporal models. We extend this idea to temporal models, where we assume that each feature is represented by a time series of length L and that the time series are aligned in time, so that the data is represented by an M × D × L tensor, with instance i represented by a matrix X (i) and feature j of instance i by a time series x By unrolling in time, this may be viewed as tabular data consisting of M instances and D·L features, so that permutations of individual features in the tabular setting correspond to permutations of individual timesteps in the temporal setting. However, doing so ignores the temporal structure of the data and correlations within time series. Thus, we consider joint permutations of timesteps in the temporal setting.
While joint permutations of subsets of features have been studied for tabular data (Henelius et al., 2014;Gregorutti et al., 2015), the exponentially large number of subsets make the task intractable in general, and high-value subsets must be selected, such as through domain knowledge (Lee et al., 2019). Several removal-based explanation methods may be viewed as operating on high-value feature subsets (Covert et al., 2020a). In the temporal setting, the structure of the data lends itself to high-value subsets constituted by contiguous regions, i.e., windows, in time. The perturbed output of the model for instance i w.r.t. feature j and time window [k 1 , k 2 ] is given by: [k1,k2] is the time series for instance i and feature j with timesteps in the window [k 1 , k 2 ] replaced by the corresponding window from another instance l = i, as shown in Figure 2b.
(7) We compute the perturbed loss L y (i) , f X (i,l) j, [k1,k2] , the change in loss (Equation 3), and average change in loss (Equation 4) for instance i. We average this over all instances i = 1 . . . M to compute the importance score corresponding to the window [k 1 , k 2 ] for feature j: By setting k 1 = 1 and k 2 = L, we permute the entire sequence of feature j and get its overall importance I (f, j, [1, L]). Equation 8.

Identifying Important Windows
Given that the features have an explicit temporal structure, we want to localize the timesteps that the model may be focusing on. We assume that for a given feature j, there exists an underlying contiguous time window W * = [k 1 , k 2 ] : 1 ≤ k 1 < k 2 ≤ L, so that most of the effect of perturbing j lies in W * . Specifically, we consider a partitioning of the sequence into three windows: prior window where W P are W S both have low importance and a size of zero or more timesteps. In order to pin down the most salient timesteps, we want to find the smallest window W * such that where the window localization parameter γ : 0 < γ < 1 controls the degree to which the model focuses on the window W * .
We use a binary search algorithm to identify the largest possible prior window W P and subsequent window W S , and by exclusion, identify the important window W * . We start with an initial estimateŴ P = [1,k 1 ] withk 1 = L 2 . We then perturbŴ P and observe its importance score I f, j,Ŵ P .
IfŴ P contains important timesteps, its importance score is likely to be inflated due to the breakage of correlations between all timesteps of the important window, i.e., predictors strongly associated with the response (Nicodemus et al., 2010), leading the search algorithm to contractŴ P to exclude these timesteps. On the other hand, ifŴ P has a low importance score I f, j,Ŵ P < 1−γ 2 I (f, j, [1, L]), we expand it unless doing so would violate this condition.
We follow a similar procedure to identify k 2 , starting from an initial estimateŴ S = [k 2 , L] withk 2 = k 1 + 1, measuring its importance score, and iteratively expanding or contracting it under the constraintk 2 > k 1 , until we identify the largestŴ S that satisfies I f, j, . We select the final boundary estimates k 1 and k 2 = L − |Ŵ S | to characterize the important window W * and compute its importance score.

Determining the Importance of Feature Ordering
Permutations of time series ordering have been used to detect circadian patterns in gene expression data (Storch et al., 2002;Ptitsyn et al., 2006). To determine the importance of the ordering of a feature j within a given window [k 1 , k 2 ], we permute its values within the window, as illustrated in Figure 2c, and average across instances. Let Π [k1,k2] be a permutation over timesteps within the window: Π [k1,k2] = Permute [ k 1 , k 1 + 1, . . . , k 2 ] = π k1 , π k1+1 , . . . , π k2 . The perturbed model output is given by: where the permuted time series for instance i and feature j is given by: We compute the change in loss between the perturbed and baseline loss (Equation 3) for multiple permutations Π [k1,k2] , and use hypothesis testing to determine the importance of the ordering of j.

Hypothesis Testing and False Discovery Rate Control
Several works have leveraged hypothesis testing in conjunction with permutations to ascertain feature importance (Ojala & Garriga, 2010;Mentch & Hooker, 2016;Lee et al., 2019;Burns et al., 2020). We use hypothesis testing to evaluate feature importance, while using the scores from Equation 5 to quantify the degree of importance. Intuitively, if a feature is useful to the model for predicting the target, then permuting its values across instances should increase the model's loss on average.
Specifically, we use the formulation of permutation tests in Ojala & Garriga (2010) to evaluate the effect of permuting a feature and ascertain its importance. We use mean loss as the test statistic, so that the empirical p-value for testing the significance of permuting feature j is given by: where P j is a set of permutations of the original data with feature j permuted in some way,L is the mean loss for the original data, andL Π is the mean loss for permuted data Π.
Depending on the choice of P j , we can use Equation 11 to test the overall importance, window importance, and ordering importance for a feature j. These tests may be organized as a hierarchy, as shown in Figure 3a, so that a test is performed only if its parent test, if any, returns a significant p-value.
The multiplicity of hypothesis tests for a given feature and across the set of features leads to a multiple comparisons problem. We address this by using a hierarchical false discovery rate (FDR) control methodology (Yekutieli, 2008), with the FDR for tests with the same parent controlled using the Benjamini-Hochberg procedure (1995). This approach also readily extends to the case when features may be arranged in a hierarchy in order to interpret models in terms of feature groups (Henelius et al., 2014;Gregorutti et al., 2015;Lee et al., 2019), as shown in Figure 3b.

Results
We evaluate TIME in two ways: by analyzing synthetic data sets and models where the ground truth pertaining to relevant features and their temporal properties is known, and by analyzing a long short term memory (LSTM) model trained to predict in-hospital mortality from intensive care unit (ICU) data.

Simulations
To evaluate our approach, we create synthetic time-series data where we control the generating processes for different features and the functions used to generate targets for the instances. We also create synthetic models that approximate these functions and serve as the models to be analyzed.
We control the features that are relevant to the models, as well as the temporal properties of the models, including the relevant window and dependence on ordering for each feature. We then analyze these models using TIME and evaluate the results in terms of power (the fraction of rele-vant features correctly identified) and FDR (the fraction of features estimated to be important, but not truly relevant in the underlying function).
Synthetic data and models. We use random Markov chains to generate time series data with both categorical and continuous features, and randomly selected windows for each feature. More details of the synthetic data generation model are given in Section S2. For each synthetic data set, we create a multi-level function to generate targets for the instances. For each feature j, we apply a feature function g j that aggregates the values within the window [k 1 , k 2 ] of that feature: where (i)ḡ j is aggregation operator, randomly selected from one of max, average (both insensitive to temporal ordering), monotonic-weighted-average and random-weighted-average (both sensitive to temporal ordering) functions, (ii)g j is randomly selected from one of identity, absolute-value and square functions and serves to potentially induce interactions between the timesteps in the window, and (iii)ĝ j is a standardization operator that yields zero mean and unit variance for feature j across the data set.
We designate a subset of all features as relevant and take a linear combination of their feature functions to generate the target: where R is the set of relevant features, and coefficients α j are sampled uniformly at random between -1 and 1. This serves to generate responses for a regression task. To emulate a classification task, we choose a threshold such that half the instances are labeled negative and the other half are labeled positive.
The synthetic model represents an approximation of this function and is generated by adding a weighted linear combination of the set of irrelevant features R to the function: The terms corresponding to the irrelevant features represent noise in the model, with the overall level of noise controlled by the multiplier β. The rationale behind the approximation is to have a realistic model that does not perfectly match the underlying function and whose output changes in a small way when irrelevant features are perturbed, but not in a way that consistently affects the loss function L.
Unlike a real model where training may involve optimizing over a loss function, here we use a loss function only to measure the fidelity of the model output f to the target y and compute the importance of each feature. We use quadratic loss for regression models and binary cross-entropy for classification models.
Baseline comparisons. We compare TIME against three baseline methods: LIME (Ribeiro et al., 2016): a model-agnostic method for local explanations. We aggregate local feature importance scores to generate global ones, based on the submodular pick algorithm described by the authors. We include LIME due to its widespread usage as an explanation method, and as a representative of other methods that focus on the model output rather than loss and generate local explanations.
SAGE (Covert et al., 2020b): a model-agnostic method that generalizes SHAP (Lundberg & Lee, 2017) to global explanations. SAGE is intractable to compute exactly, so we use two approximations provided by the authors: sampling held-out features from (i) their marginal distributions, or (ii) from reference values (for which we use mean values); we refer to these as SAGE and SAGE-m respectively.
PERM: a method using traditional permutation tests on individual, i.e., tabular features, rather than sequences. This is also a global feature importance method, although it does not take into account the temporal structure of the data.
Since the baseline methods are designed for a tabular feature representation, we unroll the temporal data comprising D features and L timesteps into tabular data with D · L features. For the sake of brevity, and to avoid confusion with temporal features, we refer to tabular features simply as 'timesteps' in the context of evaluation, since each tabular feature corresponds to a single feature-timestep pair in the original representation.
For TIME, we set the window localization parameter γ to 0.99 and control the false discovery rate at 0.1. We sample 50 permutations to compute each p-value.
We generate data sets with 1,000 instances, 10 features and 20 timesteps per feature. Five features are randomly selected as relevant. We create a synthetic model for each data set, with β tuned to yield a 90% accuracy (for classification models) or an R 2 value of 0.9 (for regression models). We evaluate the methods by examining power and FDR for identifying relevant features as well as timesteps, and average the results over 100 data sets and models.
For the baseline methods, we estimate a feature's importance by averaging non-zero importance scores across the timesteps belonging to feature. We sort timesteps in decreasing order of importance scores and report the n features or We refer to these as TIME and TIME-n respectively. Table 1 shows results from this comparison, averaged across 100 data sets and classification models. Both TIME and TIME-n provide higher power and lower FDR than all baselines for both features and timesteps, and the average FDR is well-controlled at the 0.1 level. Figure 4 illustrates feature importance explanations for one such model. It shows a set of heat maps indicating relevant timesteps for the ground truth model along with the importance scores returned by TIME and SAGE. For the ground truth model, boxes corresponding to relevant timesteps are shown in a uniform color. For the explanation methods, colored boxes indicate non-zero importance scores, with higher scores shown in darker shades. Hatched textures are used to show features for which ordering is relevant (ground truth) or identified as important (TIME), but they are not shown for SAGE since it is not able to detect the significance of ordering. TIME assigns importance scores to windows for each feature, while SAGE (as well as other baseline methods) assign importance scores to each timestep, since they operate on a tabular representation. For this model, TIME identifies all the relevant features, timesteps and their ordering correctly. SAGE assigns non-zero importance scores to all the relevant timesteps, but in some cases, irrelevant timesteps are ranked above relevant ones, adversely affecting its power and FDR for detecting important features.
Performance vs. test set size. In addition to baseline comparisons, we examine the performance of our method as a function of the size of the test set used to analyze the model. We generate data sets with 6,400 instances, 30 features and 50 timesteps per feature, and increase the size of the test set available to the model in multiples of two. Ten features are randomly selected as relevant. For each test set size, we aggregate the results over 100 different models. Figure 5 shows the results of this analysis for regression models, and similar results are obtained for classification models. Figure 5a shows average power and FDR for relevant features and timesteps as a function of test set size. The power increases as the test set size increases and has high terminal values, indicating that our approach is successful at identifying most of the relevant features and windows. The average FDRs are well-controlled at the 0.1 level. Figure 5b shows average power and FDR for detecting features and windows for which the ordering of values is im-portant. Feature ordering refers to the ordering of a feature's values across its entire sequence. Since the distribution of values inside the window is different from that outside the window, the model is sensitive to the ordering of all features having windows smaller than the sequence length. However, the model is sensitive to the ordering of values within the window only for certain feature functions. At the largest test set size TIME is able to detect ordering with high accuracy while FDRs are well-controlled at the 0.1 level. We detect window ordering at lower power compared to feature ordering due to the greater difficulty of the task, and the fact that relevant features that are not identified as important are not assessed for important windows or their ordering.

MIMIC-III Benchmark LSTM Model
To consider a challenging, real-world task, we analyze an LSTM model trained on the MIMIC-III database, a publicly available critical care database consisting of records of 58,976 intensive care unit (ICU) admissions (Johnson et al., 2016). The model is one of several proposed as part of a benchmark suite for four different clinical prediction tasks over the MIMIC-III database (Harutyunyan et al., 2017). In this analysis, we examine a regular LSTM (Hochreiter & Schmidhuber, 1997) trained to predict in-hospital mortality of patients given the first 48 hours of their ICU stay observations. After pruning to meet inclusion criteria (such as stays of 48 hours or longer), the data comprises training, validation and test sets of 14,682, 3,221 and 3,236 stays respectively, with 13.23% of the labels being positive. There are 76 features, each represented by a sequence of length 48. The features are derived from chart and laboratory measurements, and include 'mask' features indicating interpolated values. Further details on the model and features may be found in the benchmarking paper (Harutyunyan et al., 2017).
We use the validation set to analyze the LSTM and identify important features and windows, and whether or not their ordering is important to the model. We set the window localization parameter γ as 0.95 and control the false discovery rate at 0.1. We sample 200 permutations to compute each Each row corresponds to a single feature and shows the window corresponding to important timesteps in color. The importance score is indicated by the color bar, with higher scores associated with darker shades, and hatched textures show windows that were found to be significant in relation to ordering.
p-value. Figure 6 shows the results of this analysis. The feature importance analysis of the MIMIC-III model identifies a set of 31 features that are important for the model's predictions, as well as the important windows for these features. The important windows almost always focus on the more recent part of the patients' histories, which is expected since death is more likely to be indicated by abnormalities in the later stages of the ICU stay. We also note that the ordering of timesteps is found to be important for some features, suggesting that the model may be picking up on trends for these features.
To validate our analysis, we use the set of features and windows estimated to be important to perform feature selection. We prune the features that are not estimated to be important and set the out-of-window timesteps for important features to zero. We then retrain the LSTM on the pruned data set and compare its performance to the original model on the held-aside test set. We further train and test 20 featureselected models with 31 features and windows chosen at random. The area under the ROC curve for the retrained model on the test set is 0.836, which is close to that of the original model (0.838) but significantly higher than the models using randomly selected features (mean 0.801, SD 0.015), suggesting that TIME is able to identify a salient subset of features and windows for this model.

Conclusions
We have presented TIME, a method to explain black-box models having an explicit sequential or temporal structure. TIME identifies the set of important features and their de-gree of importance, and for each important feature, it identifies the window that the model focuses on and the significance of its ordering. It uses hypothesis testing and an FDR control methodology to detect these with statistical rigor.
Our experiments show that on synthetic data, TIME performs significantly better than baseline methods at identifying relevant features and timesteps. Additionally, in comparison to baselines, TIME is efficient to compute and potentially more interpretable since it identifies contiguous windows rather than scattered timesteps. We apply TIME to an LSTM trained to predict risk of in-hospital mortality from ICU data, and we identify salient features, windows and ordering in patients' clinical histories that the model focuses on. We show that a model trained using features and timesteps selected using this analysis performs nearly as well as the original model.
There are several limitations of this work that we plan to address in future research. Our approach for permutation across sequences currently assumes regularly sampled, timealigned and fixed-length sequences, which excludes models trained over irregularly sampled, variable-length sequences. We assume that there exists a single contiguous window that is important, which may not be the case. We plan to explore the use of conditional permutation tests (Strobl et al., 2008) to avoid breaking across-feature correlations, though efficient and accurate estimation of conditional distributions for temporal data remains a challenging task. We also plan to extend our approach to analyze sequence-sequence models, and to consider windows that are aligned in other ways, such as on an absolute scale (e.g., dates on the Gregorian calendar) or a relative scale (e.g., patient age).

S1. Importance Scores
Equation 8, reiterated below, computes the importance score for window [k 1 , k 2 ] of feature j: Intuitively, this captures the notion that a window for a given feature is considered important to the model if it has a positive association with the target, in the sense that permuting its values increases the model's loss on average. In Proposition S1.1, we show that for additive models and under certain assumptions, Equation 8 captures this association in terms of the covariance between the target and the feature function operating on the window. Proposition S1.1. Let X 1 . . . X D be independent random vectors of size L representing temporal features F = {1, 2, . . . , D} for the additive model f (X 1 , . . . X D ) = j ∈F g j (X j ). Let W be a window for X 1 perturbed according to Equation 7, and let W represent the timesteps outside the window. Let g 1 (X 1 ) decompose additively over the sequence as: where X j,W and X j,W represent subsequences of X 1 inside and outside the window respectively, and g j,W (X j,W ) and g j,W (X j,W ) represent functions over these subsequences. Let Y be the target and L be the quadratic loss function.
Then, the importance score I(f, X 1 , W ) of the window W for feature X 1 as computed by Equation 8 satisfies: E [I(f, X 1 , W )] = 2 cov Y, g j,W (X j,W ) − cov g j,W (X j,W ), g j,W (X j,W ) . (S1) Proof. Our formulation of the importance score corresponds to an empirical estimate of difference-based model reliance (Fisher et al., 2019), and equivalently, when quadratic loss is assumed, the permutation importance score in Gregorutti et al. (2017), extended to the temporal setting.
In the tabular setting, when the model is assumed to be additive, the expected overall importance score of feature X 1 is given by (Fisher et al., 2019): where Xj corresponds to the subset of features other than j. Extending this to temporal models, the expected importance score for window W of feature X 1 is given by: since cov g j,W (X j,W ), g j (Xj) = 0 ∀j ∈ F\j as X 1 , . . . X D are independent.

S2. Synthetic Data Generation
We use Markov chains to generate time series data, as shown in Figure S1a. Each feature is associated with a randomly selected window and a pair of Markov chains, one each to generate values for in-window and out-of-window timesteps. The number of states in each chain is sampled uniformly at random between 2 and 5. The features include a combination of continuous and discrete features. Each state m is associated with a Gaussian random variable S m ∼ N (µ m , σ 2 m ) (for continuous features) or an integer value (for discrete features) and transition probabilities p mn to other states n within the same chain. The mean, standard deviation, and transition probabilities for each state are sampled uniformly at random. The sequence for a given instance and feature is generated via a random walk through the chains. For example, a sequence i for feature j with 5 timesteps and underlying window [2, 3] may be generated as: j,5 = s 0,1 , s 2,2 , s 1,3 , s 1,4 , s 0,5 , where s m,t is sampled from S m at timestep t. For some continuous features, sampled values are aggregated over time to model increasing, constant, or decreasing trends, as shown in Figure S1b. In this case, for each timestep t, x j,t = s m,t + t−1 k=1 x j,k . Feature Importance Explanations for Temporal Black-Box Models  Figure S1. (a) Generator for a continuous feature consisting of two Markov chains, one each for in-window and out-of-window states. Here each Markov chain consists of three states, and each state is associated with a Gaussian random variable. (b) Three sequences generated via random walks through the chains, with the sampled values aggregated over time to create trends. The window is represented by blue shading.

S3. Empirical Results
Additional baseline comparisons. In addition to the comparisons shown in Table 1, we also tested variants of two baseline methods: (i) SAGE-z, an approximation of SAGE provided by its authors where held-out features are sampled from all-zero reference values, instead of mean values as in SAGE-m, and (ii) PERM-f, a variant of PERM where FDR is controlled at the 0.1 level (as is done for TIME) using the BH procedure (Benjamini & Hochberg, 1995) after computing importance scores for each timestep. These results are shown in Table S1 along with the results  from Table 1. We also performed baseline comparisons using a larger feature set composed of 30 features, 10 out of which were randomly selected as relevant. Table S2 shows these results, aggregated over 100 different models. In case of SAGE and SAGE-m, these results are aggregated over 86 and 99 models respectively, as convergence was not achieved for some models within the designated time limit (72 hours per model). PERM-f is not included since it did not identify any features as important. The results corroborate the conclusions drawn from Table 1. Performance vs. test set size for classification models. Figure S2 shows results for synthetic classification models analogous to those for synthetic regression models shown in Figure 5. The results show the same patterns as for synthetic regression models: the power for detecting relevant features, timesteps and ordering increases as the size of the test set used to examine the model is increased, and the FDR is well-controlled at the 0.1 level.

S4. Computation
Time complexity. The time complexity of our method is O(M P L min{R log L, D}), where M , D, L, R and P are the test-set size, number of features, sequence length, number of relevant features and number of permutations respectively. The logarithmic term arises since the window search algorithm partitions the search space in half at each step of the window search. In practice, as shown in Table 1, TIME is often significantly faster than PERM, since the permutation of high-value subsets (windows) can leverage a vectorized implementation to yield faster performance than permutations of their constituent timesteps.
Distributed implementation. TIME can be computed significantly faster by leveraging a distributed computing Feature Importance Explanations for Temporal Black-Box Models environment, where each node analyzes a feature or a subset of features. The results shown in Figure 5 and Figure 6 were produced with a distributed implementation using HTCondor. For a fair comparison, the distributed implementation was disabled while comparing running times of different baseline methods shown in Table 1.