You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: supplementary_examples.Rmd
+15-14Lines changed: 15 additions & 14 deletions
Original file line number
Diff line number
Diff line change
@@ -74,7 +74,7 @@ library(dplyr) # work with data frames
74
74
library(ggplot2) # plots
75
75
library(patchwork) # combine plots
76
76
77
-
# only used for heart beat detection (dependencies of `find_abp_beats()`)
77
+
# Packages only used for heart beat detection (dependencies of `find_abp_beats()`)
78
78
# install.packages("RcppRoll")
79
79
# install.packages("purrr")
80
80
@@ -86,7 +86,7 @@ source("functions.R")
86
86
87
87
# Pulse pressure
88
88
89
-
The first example corresponds to the model shown in Fig. 2 in the paper. The aim is to model variation in pulse pressure using information about the respiratory cycle.
89
+
The first example corresponds to the model shown in Fig. 2 in the paper. The aim is to model variation in pulse pressure using information about the respiratory cycle. The model is then used to estimate pulse pressure variation (PPV).
To find the extrema of the smooth, we can generate a grid of predictions using only the one smooth term `s(insp_rel_index)`. We could use `predict(type="terms")` but `gratia::smooth_estimates()` conveniently returns the values of the smooth over the original range of the independent variable (here `insp_rel_index`).
341
+
where $\alpha$ is the model constant (intercept). To find the extrema of the smooth, we can generate a grid of predictions using only the one smooth term `s(insp_rel_index)`. We could use `predict(type="terms")` but `gratia::smooth_estimates()` conveniently returns the values of the smooth over the original range of the independent variable (here `insp_rel_index`).
342
342
343
343
```{r}
344
344
insp_rel_index_smooth <- smooth_estimates(PP_gam,
@@ -355,7 +355,7 @@ sprintf("PPV is %.1f%%", PPV_est*100)
355
355
This PPV is an estimate from the model. We should also report the uncertainty. `mgcv` lets us sample parameters from the posterior distribution of the model (similarly to posterior sampling from a Bayesian model; see [Gavin Simpson's answer on StackOverflow](https://stats.stackexchange.com/questions/190348/can-i-use-bootstrapping-to-estimate-the-uncertainty-in-a-maximum-value-of-a-gam) for an in-depth example of how this can be done). Here we need posterior samples of a smooth, and, conveniently, `gratia::smooth_samples()` does just that. For each sampled smooth, we can calculate PPV, and use the many samples of PPVs to calculate a confidence interval for PPV (this approach neglects that the model intercept is also an estimate, but since the standard error for the intercept is *very* small, this has negligible effect on the width of the confidence interval).
356
356
357
357
```{r}
358
-
#| fig.cap = "Plot of 50 of the 5000 sampled smooths from the posterior distribution of the respiratory cycle smooth s(insp_rel_index)."
358
+
#| fig.cap = "Plot of 50 of the 5000 sampled smooths from the posterior distribution of the respiratory cycle smooth `s(insp_rel_index)`."
359
359
set.seed(1)
360
360
# Sample 5000 smooths
361
361
insp_smooth_samples <- smooth_samples(PP_gam, term = "s(insp_rel_index)", n = 5000)
#| fig.cap = "Histogram of PPVs calculated from each of the 5000 sampled smooths from the posterior distribution of the respiratory cycle smooth s(insp_rel_index)."
370
+
#| fig.cap = "Histogram of PPVs calculated from each of the 5000 sampled smooths from the posterior distribution of the respiratory cycle smooth `s(insp_rel_index)`."
371
371
# A function that returns PPV given a smooth and an intercept
First, we need to calculate each CVP sample's position in both the cardiac and respiratory cycles. For the respiratory cycle, we will fit a cyclic spline based on the relative position (as in the pulse pressure example above). For the cardiac cycle, we cannot simply use a cyclic spline, as the cycles vary in length. We could use a cyclic spline based on the relative position in the cardiac cycle of a CVP sample, but that would assume that the CVP waveform of a long cardiac cycle is simply a stretched version of a short cycle. Instead we assume that the cardiac cycle effect depends on the time since the P-wave (initiation of atrial contraction), without constraining it to be cyclic. Instead of trying to detect P-waves from the ECG, we assume that they appear a constant interval before the QRS complex.
445
+
First, we need to calculate each CVP sample's position in both the cardiac and respiratory cycles. For the respiratory cycle, we will fit a cyclic spline based on the relative position (as in the pulse pressure example above). For the cardiac cycle, we cannot simply use a cyclic spline, as the cycles vary in length. We could use a cyclic spline based on the relative position in the cardiac cycle of a CVP sample, but that would assume that the CVP waveform of a long cardiac cycle is simply a stretched version of a short cycle. Instead we assume that the cardiac cycle effect depends on the time since the Pwave (initiation of atrial contraction), without constraining it to be cyclic. Instead of trying to detect Pwaves from the ECG, we assume that they appear a constant interval before the QRS complex.
446
446
447
447
In the sample data, QRS complexes have already been detected from `sample_cvp$ecg`. This was done with [`rsleep::detect_rpeaks()`](https://cran.r-project.org/package=rsleep).
448
448
@@ -451,17 +451,18 @@ To find the PQ interval, we align 30 seconds of ECG recording by the detected QR
ggplot(aes(pre_qrs_index-0.3, ECG_II, group = pre_qrs_n)) +
457
+
ggplot(aes(before_qrs_index-0.3, ECG_II, group = before_qrs_n)) +
458
458
geom_line(alpha = 0.3) +
459
-
scale_x_continuous(breaks = seq(-0.2, 0.4, by = 0.05))
459
+
scale_x_continuous(breaks = seq(-0.2, 0.4, by = 0.05)) +
460
+
labs(x = "Time - QRS-time [seconds]")
460
461
```
461
462
462
-
We can see that the P-wave starts ~150 ms before the QRS complex.
463
+
We can see that the Pwave starts ~150 ms before the QRS complex.
463
464
464
-
In this example we will fit a GAM to the last 30 seconds before fluid administration starts (the green area in Figure \@ref(fig:cvp-overview)). We add each sample's position in both the cardiac cycle (starting at the P-wave) and the respiratory cycle, and filter the data to the relevant section.
465
+
In this example we will fit a GAM to the last 30 seconds before fluid administration starts (the green area in Figure \@ref(fig:cvp-overview)). We add each sample's position in both the cardiac cycle (starting at the Pwave) and the respiratory cycle, and filter the data to the relevant section.
465
466
466
467
```{r}
467
468
PQ_interval <- 0.150 #seconds
@@ -478,7 +479,7 @@ cvp_df_30 <- filter(cvp_df,
478
479
head(cvp_df_30)
479
480
```
480
481
481
-
In the above table, `P_wave_index` is the position in the cardiac cycle (seconds since latest P-wave) and `insp_rel_index` is the position in the respiratory cycle (time since inspiration start relative to length of the respiratory cycle).
482
+
In the above table, `P_wave_index` is the position in the cardiac cycle (seconds since latest Pwave) and `insp_rel_index` is the position in the respiratory cycle (time since inspiration start relative to length of the respiratory cycle).
482
483
483
484
Before fitting the model, we can visualise (individually) the three effects we subsequently want to model (collectively): time, respiratory cycle and heart cycle.
0 commit comments