Supplementary Material for: "Learning a Causal Model for Intracranial Pressure in Patients with Traumatic Brain Injury"

1 Supplementary Material A How to Read a Conditional Intensity Matrix Reading a CIM can be tricky. For example, let Gbe a directed graph with three vertices V={X, Y, Z}, connected as depicted in Figure 1. Y Z X Fig. 1: The graph G. QXi|Π Π Πi= q311 q312 . . . q31n q321 q322 q32n . . .... q3n1q3n2q3nn q211 q212 . . . q21n q221 q222 q22n . . .... q2n1q2n2q2nn q111 q112 . . . q11n q121 q122 q12n . . .... q1n1q1n2q1nn q011 q012 . . . q01n q021 q022 q02n . . .... q0n1q0n2q0nn Π Π Πi Xt+δt i Xt i Fig. 2: Generic CIM QXi|Π Π Πirepresentation. If Gis the graph of a CTBN, then each vertex in Vis mapped to a random variable in X, which in turn is parametrized by a CIM. Figure 2represents the layout of a generic CIM QXi|Π Π Πi. In this example, the three random variables have the following domains: X={0,1,2,3},Y={0,1}and Z={0,1}. When we set Y= 0 and Z= 1 the corresponding IM QX|{Y=0,Z=1}is represented as depicted in Figure 3. (0, 0) (0, 1) (1, 0) (1, 1) Conditioned 0 1 2 3 To State 0 1 2 3 From State -4.635 0.396 1.242 1.682 1.906 -2.353 1.445 0.503 1.491 0.210 -4.630 0.445 1.237 1.746 1.943 -2.631 -2.671 1.263 0.967 1.226 1.097 -2.714 1.592 0.188 0.921 0.655 -3.635 1.254 0.653 0.796 1.077 -2.668 -5.474 0.679 0.332 0.592 1.903 -3.015 1.041 1.359 1.935 1.400 -3.200 0.692 1.636 0.936 1.828 -2.643 -3.966 1.885 0.268 0.838 0.451 -4.973 0.472 0.616 1.942 1.236 -1.459 1.675 1.573 1.852 0.718 -3.129 0123 To State 0 1 2 3 From State -2.671 1.097 0.921 0.653 1.263 -2.714 0.655 0.796 0.967 1.592 -3.635 1.077 1.226 0.188 1.254 -2.668 5 4 3 2 1 0 1 Fig. 3: Specific IM QX|{Y=0,Z=1}representation. 2 If we set X= 1 then the row QX=1|{Y=0,Z=1}encodes both the exponential distribution modeling the expected transition time and the categorical distribution of the next state. Figure 4depicts in blue qiand the corresponding probability density function (PDF) of the exponential Exp(qi), while in red θij =qij/qi, i =j and the corresponding distribution Cat(θ θ θi). 0123 To State 0 1 2 3 From State -2.671 1.097 0.921 0.653 1.263 -2.714 0.655 0.796 0.967 1.592 -3.635 1.077 1.226 0.188 1.254 -2.668 Diagonal ( X = 0, X = 0) 0123 Time 0.0 0.5 1.0 1.5 2.0 2.5 PDF Exponential for X = 0 = 2.671 E(X) = 0.3744 0123 To State 0 1 2 3 From State -2.671 1.097 0.921 0.653 1.263 -2.714 0.655 0.796 0.967 1.592 -3.635 1.077 1.226 0.188 1.254 -2.668 Row X = 0 0123 Category 0.0 0.2 0.4 0.6 0.8 1.0 Probability Categorical for X = 0 0123 To State 0 1 2 3 From State -2.671 1.097 0.921 0.653 1.263 -2.714 0.655 0.796 0.967 1.592 -3.635 1.077 1.226 0.188 1.254 -2.668 Diagonal ( X = 1, X = 1) 0123 Time 0.0 0.5 1.0 1.5 2.0 2.5 PDF Exponential for X = 1 = 2.714 E(X) = 0.3685 0123 To State 0 1 2 3 From State -2.671 1.097 0.921 0.653 1.263 -2.714 0.655 0.796 0.967 1.592 -3.635 1.077 1.226 0.188 1.254 -2.668 Row X = 1 0123 Category 0.0 0.2 0.4 0.6 0.8 1.0 Probability Categorical for X = 1 0123 To State 0 1 2 3 From State -2.671 1.097 0.921 0.653 1.263 -2.714 0.655 0.796 0.967 1.592 -3.635 1.077 1.226 0.188 1.254 -2.668 Diagonal ( X = 2, X = 2) 012 Time 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 PDF Exponential for X = 2 = 3.635 E(X) = 0.2751 0123 To State 0 1 2 3 From State -2.671 1.097 0.921 0.653 1.263 -2.714 0.655 0.796 0.967 1.592 -3.635 1.077 1.226 0.188 1.254 -2.668 Row X = 2 0123 Category 0.0 0.2 0.4 0.6 0.8 1.0 Probability Categorical for X = 2 0123 To State 0 1 2 3 From State -2.671 1.097 0.921 0.653 1.263 -2.714 0.655 0.796 0.967 1.592 -3.635 1.077 1.226 0.188 1.254 -2.668 Diagonal ( X = 3, X = 3) 0123 Time 0.0 0.5 1.0 1.5 2.0 2.5 PDF Exponential for X = 3 = 2.668 E(X) = 0.3748 0123 To State 0 1 2 3 From State -2.671 1.097 0.921 0.653 1.263 -2.714 0.655 0.796 0.967 1.592 -3.635 1.077 1.226 0.188 1.254 -2.668 Row X = 3 0123 Category 0.0 0.2 0.4 0.6 0.8 1.0 Probability Categorical for X = 3 Title Fig. 4: Exponential and categorical distributions for QX=1|{Y=0,Z=1}. For example, here qii =−2.714 means that the parameter of the exponential distribution is λ=qi=−qii = 2.714, with an expected transition time of E[Exp(λ)] = 1/λ = 0.3685. Let’s say time is measured in hours, then with {Y= 0, Z = 1}the variable X= 1 is expected to transition out of the current x1in 0.3685 hours, or 22.11 minutes. Moreover, Xwill transition to 0with probability qij /qi= 1.263/2.714 = 0.4653. 3 B Parameter Learning from Incomplete Trajectories In real world settings it is unlikely to fully observe every variable all the time. Expectation-Maximization. A more realistic approach is to estimate the expected sufficient statistics Nijk and Tij from an incomplete trajectory σ−via the expectation-maximization (EM) [1] algorithm: 0. Initialize the current parameters Q0 Xi|Π Π Πiarbitrarily. 1. Expectation step: with the current k-th parameters, use σ−as evidence to generate a completed trajectory σ+and compute Nijk and Tij from σ+. 2. Maximization step: estimate the next (k+ 1)-th parameters from Nijk and Tij as if they were computed from a complete trajectory. 3. Check convergence using a stopping criterion, e.g. absolute difference from k-th and (k+1)-th parameters or maximum number of iterations. If no convergence is reached, then go to (1), otherwise return the current parameters. The difficult part of the EM algorithm is the expectation step: completing a partially-observed trajectory is not a trivial task. Importance Sampling. A method to sample a completed trajectory σ+from the evidence econtained in the incomplete trajectory σ−is given by the Importance Sampling [2] algorithm, which generates a set of weighted trajectories compatible with e, as depicted in Figure 5. Since the trajectories are drawn from a proposal distribution P′consistent with e, each σ+ iis weighted by the likelihood of the evidence ω(σ+ i) = P(σ+ i,e)/P′(σ+ i), and the expected sufficient statistics are computed as a weighted average over ω(σ+ i). Point Trans. Interval σ − σ+ 1 σ+ 2 x0 x1 x0 x1 x0 x1 0 1 2 3 4 5 Fig. 5: Evidence coming from the incomplete trajectory σ−can be point-wise, transition-like or interval-like. More than one completed trajectory σ+generated via importance sampling, such as σ+ 1and σ+ 2, is compatible with the same σ−. 4 Uncertain Evidence. In the context of event-based data, evidence comes is many forms. I can happen that there is uncertainty about the exact status in which a variable is at a specific point in time due contrasting information coming from different data sources. To deal with this kind on uncertainty even the evidence must be treated as uncertain. Luckily, importance sampling is flexible enough to accommodate for uncertain evidence [3] without changing the whole theoretical framework. C Structure Learning from Incomplete Trajectories As EM learn the parameters of a CTBN from incomplete trajectories, the only remaining piece to construct a CTBN from data is to learn its graph. Structural Expectation-Maximization. Learning a graph from incomplete trajectories can be done via the structural EM (SEM) [4] algorithm: 0. Initialize the current graph G0and parameters Q0 Xi|Π Π Πiarbitrarily. 1. Expectation step: with the current k-th graph and parameters, use σ−as evidence for the EM algorithm to generate a completed trajectory σ+. 2. Maximization step: estimate the next (k+ 1)-th graph using a σ+and the prior knowledge Kas input for a structure learning algorithm. 3. Check convergence using a stopping criterion, e.g. absolute difference from k-th and (k+ 1)-th graph and/or parameters, or maximum number of iterations. If no convergence is reached, then go to (1), otherwise return the current solution. Causal Discovery. The structure learning task can be constrained by experts’ knowledge for two distinct purposes: (i) to speed up the learning step, by discarding irrelevant edges or forcing them in place, thus skipping the independence checks altogether, and (ii) to add a causal interpretation of it, making it not only a predictive model but also a causal model, capable of explaining the interplay between the observed variables. In this sense, while structure learning aims to construct a model for its inference capability, the causal discovery task serves as a tool to investigate an understudied problem, where the resulting model make the cause-effect relationships explicit. 5 D CENTER-TBI Categorical Levels Table 1: Categorical levels of selected variables. Variable Levels HVICP –(0,15] - ICP lower or equal to 15 mmHg. –(15,20] - ICP between 15 and 20 mmHg. –(20,25] - ICP between 20 and 25 mmHg. –(25,+∞)- ICP higher than 15 mmHg. HosComplEventHypocapnia –0 - No episodes, –1 - Single episode, –2 - Multiple episodes. MarshallCTClassification –1 - No visible pathology on CT, –2 - Cisterns present, MLS < 5 mm, –3 - Cisterns compressed or absent, MLS < 5 mm, –4 - MLS > 5 mm, no mass lesion > 25 cc, –5 - Evacuated mass lesion, –6 - Non-evacuated mass lesion. TILCSFDrainage –0 - No, –1 - Yes. TILFever –0 - No, –1 - Yes. TILHyperosmolarTherapy –0 - No, –1 - Yes. TILHyperventilation –0 - No, –1 - Yes. TILMannitolDose –0No mannitol dose. –(0,2] - Total mannitol dose lower or equal to 2g. –(2,+∞)- Total mannitol dose higher than 2g. TILSedation –0 - No, –1 - Yes. 6 E Prior Knowledge Elicitation To guide the causal discovery procedure, prior knowledge on the interaction mechanism between ICP and treatments is encoded as required edges in Table 2. Table 2: Required edges elicited from experts knowledge. Required Edges Description HVICP ⇆ HosComplEventHypocapnia Elevated ICP is often managed with hyperventilation, which lowers PaCO2 and can lead to hypocapnia; hypocapnia causes vasoconstriction, reducing ICP. HVICP ⇆ MarshallCTClassification Severe CT abnormalities are associated with higher ICP; persistent intracranial hypertension can also worsen injury progression detectable on CT. HVICP ⇆ TILCSFDrainage CSF drainage is used to lower ICP by reducing intracranial volume; in turn, sustained ICP elevation is a trigger for CSF drainage. HVICP ⇆ TILFever Fever raises cerebral metabolism and blood flow, which can increase ICP; uncontrolled ICP prompt aggressive management to limit metabolic demand. HVICP ⇆ TILHyperosmolarTherapy Hyperosmolar therapy lowers ICP by drawing water out of brain tissue; intracranial hypertension is the main indication for escalating therapy. HVICP ⇆ TILHyperventilation Hyperventilation reduces ICP via cerebral vasoconstriction; it is typically applied reactively in response to intracranial hypertension. HVICP ⇆ TILMannitolDose Mannitol is administered specifically to lower elevated ICP; dosing is driven by the severity and persistence of ICP elevation. HVICP ⇆ TILSedation Sedation reduces cerebral metabolic rate and prevents ICP surges; worsening ICP is a frequent reason to escalate sedation intensity. 7 F Informative Intensity Matrices Figure 6is associated with a patient with a severe TBI (Marshall CT 5), complicated by hypocapnia, treated with mannitol at low doses. The ICP is markedly unstable, with transitions out of the <15 mmHg range to the 20-25 mmHg or ≥25 mmHg ranges. This suggests that with hypocapnia present, the therapeutic regime is insufficient to avoid ICP excursions. <15 15-20 20-25 >=25 To State <15 15-20 20-25 >=25 From State -1016.358 0.005 1016.347 0.005 66.231 -66.234 0.002 0.002 0.001 0.002 -0.005 0.001 0.017 801.970 0.017 -802.003 Diagonal ( X = <15, X = <15) 0.000 0.002 0.004 0.006 0.008 0.010 Time 0 200 400 600 800 1000 PDF Exponential for X = <15 = 1016 E(X) = 0.0009839 <15 15-20 20-25 >=25 To State <15 15-20 20-25 >=25 From State -1016.358 0.005 1016.347 0.005 66.231 -66.234 0.002 0.002 0.001 0.002 -0.005 0.001 0.017 801.970 0.017 -802.003 Row X = <15 <15 15-20 20-25 >=25 Category 0.0 0.2 0.4 0.6 0.8 1.0 Probability Categorical for X = <15 <15 15-20 20-25 >=25 To State <15 15-20 20-25 >=25 From State -1016.358 0.005 1016.347 0.005 66.231 -66.234 0.002 0.002 0.001 0.002 -0.005 0.001 0.017 801.970 0.017 -802.003 Diagonal ( X = 15-20, X = 15-20) 0.000 0.025 0.050 0.075 0.100 0.125 0.150 Time 0 10 20 30 40 50 60 PDF Exponential for X = 15-20 = 66.23 E(X) = 0.0151 <15 15-20 20-25 >=25 To State <15 15-20 20-25 >=25 From State -1016.358 0.005 1016.347 0.005 66.231 -66.234 0.002 0.002 0.001 0.002 -0.005 0.001 0.017 801.970 0.017 -802.003 Row X = 15-20 <15 15-20 20-25 >=25 Category 0.0 0.2 0.4 0.6 0.8 1.0 Probability Categorical for X = 15-20 <15 15-20 20-25 >=25 To State <15 15-20 20-25 >=25 From State -1016.358 0.005 1016.347 0.005 66.231 -66.234 0.002 0.002 0.001 0.002 -0.005 0.001 0.017 801.970 0.017 -802.003 Diagonal ( X = 20-25, X = 20-25) 0 500 1000 1500 2000 Time 0.000 0.001 0.002 0.003 0.004 0.005 PDF Exponential for X = 20-25 = 0.004816 E(X) = 207.6 <15 15-20 20-25 >=25 To State <15 15-20 20-25 >=25 From State -1016.358 0.005 1016.347 0.005 66.231 -66.234 0.002 0.002 0.001 0.002 -0.005 0.001 0.017 801.970 0.017 -802.003 Row X = 20-25 <15 15-20 20-25 >=25 Category 0.0 0.2 0.4 0.6 0.8 1.0 Probability Categorical for X = 20-25 <15 15-20 20-25 >=25 To State <15 15-20 20-25 >=25 From State -1016.358 0.005 1016.347 0.005 66.231 -66.234 0.002 0.002 0.001 0.002 -0.005 0.001 0.017 801.970 0.017 -802.003 Diagonal ( X = >=25, X = >=25) 0.000 0.002 0.004 0.006 0.008 0.010 0.012 Time 0 100 200 300 400 500 600 700 800 PDF Exponential for X = >=25 = 802 E(X) = 0.001247 <15 15-20 20-25 >=25 To State <15 15-20 20-25 >=25 From State -1016.358 0.005 1016.347 0.005 66.231 -66.234 0.002 0.002 0.001 0.002 -0.005 0.001 0.017 801.970 0.017 -802.003 Row X = >=25 <15 15-20 20-25 >=25 Category 0.0 0.2 0.4 0.6 0.8 1.0 Probability Categorical for X = >=25 (('HosComplEventHypocapnia', '1'), ('MarshallCTClassification', '5'), ('TILCSFDrainage', '0'), ('TILFever', '0')) (('TILHyperosmolarThearpy', '0'), ('TILHyperventilation', '0'), ('TILMannitolDose', '0-2'), ('TILSedation', '0')) Fig. 6: ICP parameters for a patient with a Marshall CT classification of 5, one hypocapnia event and daily cumulative mannitol dose between 0-2 g. 8 Figure 7is similar to the previous one, but in this case the patient does not experience any hypocapnia. Still, rapid excursions above <15 mmHg and 15-20 mmHg are present and the sojourn time in these states is short. <15 15-20 20-25 >=25 To State <15 15-20 20-25 >=25 From State -759.555 0.000 0.000 759.555 754.906 -755.067 0.000 0.161 0.000 918.104 -918.105 0.000 0.000 0.000 742.603 -742.603 Diagonal ( X = <15, X = <15) 0.0000 0.0025 0.0050 0.0075 0.0100 0.0125 Time 0 100 200 300 400 500 600 700 PDF Exponential for X = <15 = 759.6 E(X) = 0.001317 <15 15-20 20-25 >=25 To State <15 15-20 20-25 >=25 From State -759.555 0.000 0.000 759.555 754.906 -755.067 0.000 0.161 0.000 918.104 -918.105 0.000 0.000 0.000 742.603 -742.603 Row X = <15 <15 15-20 20-25 >=25 Category 0.0 0.2 0.4 0.6 0.8 1.0 Probability Categorical for X = <15 <15 15-20 20-25 >=25 To State <15 15-20 20-25 >=25 From State -759.555 0.000 0.000 759.555 754.906 -755.067 0.000 0.161 0.000 918.104 -918.105 0.000 0.000 0.000 742.603 -742.603 Diagonal ( X = 15-20, X = 15-20) 0.0000 0.0025 0.0050 0.0075 0.0100 0.0125 Time 0 100 200 300 400 500 600 700 PDF Exponential for X = 15-20 = 755.1 E(X) = 0.001324 <15 15-20 20-25 >=25 To State <15 15-20 20-25 >=25 From State -759.555 0.000 0.000 759.555 754.906 -755.067 0.000 0.161 0.000 918.104 -918.105 0.000 0.000 0.000 742.603 -742.603 Row X = 15-20 <15 15-20 20-25 >=25 Category 0.0 0.2 0.4 0.6 0.8 1.0 Probability Categorical for X = 15-20 <15 15-20 20-25 >=25 To State <15 15-20 20-25 >=25 From State -759.555 0.000 0.000 759.555 754.906 -755.067 0.000 0.161 0.000 918.104 -918.105 0.000 0.000 0.000 742.603 -742.603 Diagonal ( X = 20-25, X = 20-25) 0.000 0.002 0.004 0.006 0.008 0.010 Time 0 200 400 600 800 PDF Exponential for X = 20-25 = 918.1 E(X) = 0.001089 <15 15-20 20-25 >=25 To State <15 15-20 20-25 >=25 From State -759.555 0.000 0.000 759.555 754.906 -755.067 0.000 0.161 0.000 918.104 -918.105 0.000 0.000 0.000 742.603 -742.603 Row X = 20-25 <15 15-20 20-25 >=25 Category 0.0 0.2 0.4 0.6 0.8 1.0 Probability Categorical for X = 20-25 <15 15-20 20-25 >=25 To State <15 15-20 20-25 >=25 From State -759.555 0.000 0.000 759.555 754.906 -755.067 0.000 0.161 0.000 918.104 -918.105 0.000 0.000 0.000 742.603 -742.603 Diagonal ( X = >=25, X = >=25) 0.0000 0.0025 0.0050 0.0075 0.0100 0.0125 Time 0 100 200 300 400 500 600 700 PDF Exponential for X = >=25 = 742.6 E(X) = 0.001347 <15 15-20 20-25 >=25 To State <15 15-20 20-25 >=25 From State -759.555 0.000 0.000 759.555 754.906 -755.067 0.000 0.161 0.000 918.104 -918.105 0.000 0.000 0.000 742.603 -742.603 Row X = >=25 <15 15-20 20-25 >=25 Category 0.0 0.2 0.4 0.6 0.8 1.0 Probability Categorical for X = >=25 (('HosComplEventHypocapnia', '0'), ('MarshallCTClassification', '5'), ('TILCSFDrainage', '0'), ('TILFever', '0')) (('TILHyperosmolarThearpy', '0'), ('TILHyperventilation', '0'), ('TILMannitolDose', '0-2'), ('TILSedation', '0')) Fig. 7: ICP parameters for a patient with a Marshall CT classification of 5 and daily cumulative mannitol dose between 0-2 g. 9 G Limitations While the resulting model takes into account several aspects and issues related to the interaction between treatments and ICP under incomplete trajectories, there are still several limitations that must be addressed. Lack of Comprehensive Validation. The most relevant limitation is the lack of systematic validation, both internal and external. Internal validation must be discussed in details with clinicians to define accurately reference outcomes and targets, e.g. elevated ICP time span minimization, predicted GOSE at 6 months, and more. External validation via other data sources is currently under evaluation. Moreover, another prominent limitation for CTBNs is related to the absence of standardized evaluation procedures and readily available software. Integration of Baseline Variables. Attempts to integrate baseline variables, static patients covariates that are observed once at ICU admission, was made, see Appendix H. Still, this attempt resulted in a disconnected graph, which is typical when evidence is sparse and uncertain. A possible solution would be to adapt some approaches from existing literature on hybrid CTBNs [5,6]. Sparsity of High-Dimensional CIMs. Another limitation is related to the sparsity of the ICP CIM due to the exponential number of treatment combinations, which hinders the scalability of the model, see Appendix I. A possible solution would be to take advantage of the clinical practice to rule out inadmissible combinations, reducing the total number of parameters. Small Weights for Completed Trajectories. A side effect of sample-based inference techniques is that the mass of the completed trajectories tends to degenerate due the the sparsity of the evidence. There are two possible solutions: (i) apply mass correction [3] to the sampled trajectories to filter particles with low probability, or (ii) shift to temporally uncertain evidence that is able to model not only the uncertainty on the state, but also the uncertainty on the temporal interval, moving from a point evidence to an interval-like evidence. The former is computationally intensive, while the latter requires to manually elicit prior knowledge on the interval duration, if any.