Outline

Motivation

• Multi-center study of critical care patients from 457 ICUs ( $\approx 10k$ patients)

• maximum follow up of 60 days (we only consider short term survival $t\le 30$)

• Various confounders:

  • Age, Gender, BMI
  • Diagnosis, Admission Category
  • year of ICU admission
  • Apache II Score
  • ICU random effect

• 11-day nutrition protocol

  • prescribed calories (determined at baseline $t=0$)
  • daily caloric intake
  • daily caloric adequacy (CA) = caloric intake/prescribed calories

Caloric Intake

Motivation

• We are interested in how artificial nutrition (exposure) affects short term survival (outcome)

• Difficulty:

  • effect of nutrition might have a temporal delay (e.g. nutrition today affects survival 4 days later)

  • effect of nutrition might "wear off" after some time (e.g. nutrition on day 1 likely won't affect the hazard on day 30)

  • the (delayed) effect of nutrition also depends on the amount of nutrition (caloric adequacy) provided, possibly non-linearly

  • the same amount of exposure might have a different effect depending on the follow up and exposure time

  • the effect may be cumulative (i.e., 5 days of malnutrition in a row may be worse than only 2 in a row or 5 days malnutrition scattered throughout the follow up while on the other days "correct" amount was provided)

Terminology

We use the following terminology and notation:

• Time-to-event $t$: Time at which event times are observed

• Time of exposure $t_e$: Time at which values of the exposure are observed (must not necessarily overlap temporally with $t$, measured in the same units or be in the same domain as $t$, e.g. calendar days ( $t_e$) vs. 24h periods (days) since admission to ICU $t$)

• time-varying effects (TVE): Effects of time-constant covariates (covariates observed at the beginning of the follow-up) that can vary over time $t$

• time-dependent covariates (TDC): Covariates whose values change over time. Value changes are recorded at exposure time $t_e$ (here synonymous to exposure)

• Exposure value $z\left(t_e\right)$: The value of the TDC observed at exposure time $t_e$

• Exposure history $z$: The complete history of observed values of the exposure/TDC $z=\left(z\right(t_{e,1}),z(t_{e,2}),...,z(t_{e,Q}\left)\right)$

Terminology

A general cumulative effect/Exposure-Lag-Response Association (ELRA) can be defined as

• Partial effects $h(t,t_e,z(t_e\left)\right)$: The effect of the TDC recorded at exposure time $t_e$ with value $z\left(t_e\right)$ on the hazard at follow up time $t$ (the tri-variate function $h$ is potentially non-linear in all three dimensions)

• Cumulative effect $g(z,t)$: The total (cumulated) effect of the partial effects on the log-hazard at time $t$ given exposure history $z$

Lag-Lead-Window

The integration borders can be defined more general, such that $$g(z,t)={\int }_{t-t_{\text{lag}}-t_{\text{lead}}}^{t-t_{\text{lag}}}h(t,t_e,z(t_e\left)\right)d{t}_e$$

• Lag time $t_{\text{lag}}$: The length of the delay until the TDC recorded at exposure time $t_e$ starts to affect the hazard (often $t_{\text{lag}}=0$)
• Lead time $t_{\text{lead}}$: The duration of the effect of the TDC observed at exposure time $t_e$
• $t_{\text{lag}}$ and $t_{\text{lead}}$ define the set of exposures that contribute to the cumulative effect at time $t$ as $\left\{z\right(t_e):t_e\in [t-t_{\text{lag}}-t_{\text{lead}},t-t_{\text{lag}}\left]\right\}$
• Minimal requirement: ${\int }_{t_e:t_e\le t}$
• Special case ${\int }_0^t$ follows with $t_{\text{lag}}=0$ and $t_{\text{lead}}=t$
• Example ( $t_{\text{lag}}=4$, $t_{\text{lead}}=3$):
  • The last nutrition that will enter the cumulative effect at time $t=10$ is nutrition at $t_e\le t-t_{\text{lag}}=10-4=6$, i.e. $z(t_e=6)$
  • The earliest nutrition that will contribute to the cumulative effect at time $t=10$ is nutrition at $t_e\ge t-t_{\text{lag}}-t_{\text{lead}}=10-4-3=3$

Lag-Lead-Window

The integration borders can be defined even more general, such that

• $t_{\text{lag}}$ and $t_{\text{lead}}$ times can themselves depend on (exposure) time

• $T_e\left(t\right)$ is the set of exposure times $t_e$ relevant to the cumulative effect at time $t$

• We call $T_e\left(t\right)$ the Lag-Lead-Window or Window of effectiveness

Lag-Lead Window (Example)

ELRAs in the literature

Some models known from the literature follow as special cases of the general specification $$g(z,t)={\int }_{T_e\left(t\right)}h(t,t_e,z(t_e\left)\right)$$ when we assume that partial effects $h$ only depend on latency $t-t_e$ instead of concrete combination of $t$ and $t_e$, i.e., $h(t=30,t_e=3,z(t_e\left)\right)\stackrel{!}{=}h(t=40,t_e=13,z(t_e\left)\right)\stackrel{!}{=}\tilde{h}(t-t_e=27,z(t_e\left)\right)$

• DLNM: Distributed Lag Non-linear Models (Gasparrini et al, 2014, 2017): $g(z,t)={\int }_{T_e\left(t\right)}h(t-t_e,z(t_e\left)\right)$

• WCE: Weighted Cumulative Exposure (Sylvestre and Abrahamowicz, 2009): $g(z,t)={\int }_{T_e\left(t\right)}h(t-t_e)z\left(t_e\right)$

• Also possible within general framework:

  • more flexible WCE: $g(z,t)={\int }_{T_e\left(t\right)}h(t,t_e)z\left(t_e\right)$
  • time-varying DLNM (TV DLNM): $g(z,t)={\int }_{T_e\left(t\right)}h(t,t-t_e,z(t_e\left)\right)$

Exposure-Lag-Response Association

• $g(z,t)$ represents the cumulative, time-varying effect of exposure history $z$ on the log-hazard at time $t$

• we define its contribution to the model's additive predictor as

$$\begin{array}{c}g(z_i,t)={\int }_{T_e\left(t\right)}h({\tilde{t}}_j,t_e,z_i(t_e\left)\right)d{t}_e\approx \sum _{q:t_{e,q}\in T_e\left(t\right)}{\Delta }_{q}h({\tilde{t}}_j,t_{e,q},z_i(t_{e,q}\left)\right)\forall t\in ({\kappa }_{j-1},{\kappa }_j],\end{array}$$

with

• ${\tilde{t}}_j:=({\kappa }_j-{\kappa }_{j-1})/2,j=1,\dots ,J$
• partial effects $h({\tilde{t}}_j,t_e,z_i(t_e\left)\right)$
• quadrature weights ${\Delta }_q=t_{e,q}-t_{e,q-1}$ for numerical integration are given by the time between two consecutive exposure measurements

Tensor product smooths

Low rank representation of the tri-variate smooth function $$h(t,t_e,z(t_e\left)\right)=\sum _{\ell =1}^L\sum _{r=1}^R\sum _{m=1}^M{\gamma }_{\ell rm}{B}_m\left(z\right(t_e\left)\right)B_r\left(t_e\right)B_{\ell }\left(t\right)$$

with

• model matrix $X=X_t\odot X_{t_e}\odot X_{z\left(t_e\right)}$ and

• penalty $S={\nu }_{z\left(t_e\right)}{I}_{d_R}\otimes I_{d_L}\otimes S_{z\left(t_e\right)}+{\nu }_{t_e}{I}_{d_L}\otimes S_{t_e}\otimes I_{d_M}+{\nu }_t{S}_t\otimes I_{d_R}\otimes I_{d_M}$

$\to$ Estimate parameters $\gamma$ by optimizing $D\left(\gamma \right)+{\sum }_k{\nu }_k{\gamma }^{\prime }{S}_k\gamma$ (Wood, 2011), where

• $D\left(\gamma \right)$ is the model deviance (of the Poisson GAMM)
• $\gamma$ contains all Spline basis coefficients and random effects
• ${\nu }_k$ and $S_k,k=1,\dots ,K$ are the smoothing parameters and penalty matrices for the $k$-th smooth term, respectively

Exposure-Lag Response Association

• If we restrict the ELRA to be linear in the exposure, i.e., $h\left(z_i\right(t_e),t_e,t)=\tilde{h}(t_e,t)\cdot z_i\left(t_e\right)$ we can simplify to $$g(z_i,t)\approx \sum _{q=1}^Q{\tilde{\Delta }}_{i,q}\tilde{h}(t_{e,q},t)$$ with $${\tilde{\Delta }}_{i,q}=\{\begin{array}{cc}{z}_i\left(t_{e,q}\right){\Delta }_q& \text{if}{t}_{e,q}\in T_e\left(t\right)\\ 0& \text{else}\end{array}$$

Exposure-Lag Response Association

• Spline bases for the bivariate functions $\tilde{h}(t_e,t)$ are set up via tensor product B-spline basis with marginal bases $B_m\left(t_e\right),m=1,\dots ,M$ and $B_k\left(t\right),k=1,\dots ,K$ defined over the exposure and hazard time domains, respectively

• $M$ and $K$ delimit the maximal complexity of the ELRA

• $\tilde{h}(t_e,t)={\sum }_{m=1}^M{\sum }_{k=1}^K{\gamma }_{m,k}{B}_m\left(t_e\right)B_k\left(t\right)$

• Combining above equations yields: $$\begin{array}{c}g(z_i,t)\approx \sum _{m=1}^M\sum _{k=1}^K{\gamma }_{m,k}{\tilde{B}}_{i,m}(t_e,t)B_k\left(t\right),\end{array}$$ where ${\tilde{B}}_{i,m}(t_e,t)={\sum }_{q=1}^Q{\tilde{\Delta }}_{i,q}{B}_m\left(t_e\right)$.

Simulation - DLNM

• $\lambda \left(t|z\right)={\lambda }_0\left(t\right)\exp (\int h(t-t_e,z\left(t_e\right)\left)d{t}_e\right)$
• $t\in (0,40]$, $t_e\in [-40,40]$, $z\left(t_e\right)\in [0,10]$

Simulation (2) - TV DLNM

$\lambda \left(t|z\right)={\lambda }_0\left(t\right)\exp (\int \tilde{h}(t,t-t_e,z\left(t_e\right)\left)d{t}_e\right)={\lambda }_0\left(t\right)\exp (\int f(t)\cdot h(t-t_e,z\left(t_e\right)\left)d{t}_e\right)$ and $f\left(t\right)=-\cos \left(\pi t/t_{\text{max}}\right)$

Application

In the application example (categorical nutrition), we estimate $$\begin{array}{cc}\log \left({\lambda }_i(t|x_i,z_i,{\ell }_i)\right)& =f_0\left(t\right)+\sum _{p=1}^P{f}_p(x_{i,p},t)+g(z_i,t)+b_{{\ell }_i}\end{array}$$ with

• $f_0\left(t_j\right)={\sum }_{m=1}^M{\gamma }_{0m}{B}_m\left(t_j\right)$ represents the log baseline-hazard
• $f(x_{i,p},t_j)={\sum }_{m=1}^M{\sum }_{\ell =1}^L{\gamma }_{m\ell }{B}_m\left(x_{i,p}\right)B_{\ell }\left(t_j\right)$ are potentially non-linear, potentially non-linearly time-varying effects of confounders $x_{i,p}$
• $g(z_i,t)=g_{C2}(z_i^{C2},t)+g_{C3}(z_i^{C3},t)$
  • $z_i^{C2}$ and $z_i^{C3}$ dummy variables that indicate whether subject $i$ received category $C2$ and $C3$ nutrition on day $t_{e,q},q=1,\dots ,11$, respectively
  • $g_{C2}(z_i,t)\approx {\sum }_{q=1}^Q{\tilde{\Delta }}_{i,q}^{C2}{\tilde{h}}_{C2}(t_{e,q},t)$
• $b_{{\ell }_i}$ is the random effect associated with ICU (cluster) ${\ell }_i$ at which subject $i$ is treated

$\to C1$ reference category

PAMM

• We know how the estimate such models in the framework of
  Generalized Additive Mixed Models (GAMMs)

• Fortunately, we can fit survival models via Poisson GLMs/GAMMs by representing them as a Piece-wise exponential Additive Mixed Model (PAMMs)

• to do so requires to

  • divide the follow up $(0,t_{max}]$ into $J$ intervals with $J+1$ cut-points $0={\kappa }_0<\dots <{\kappa }_J=t_{max}$
  • transform the data into appropriate format (pseudo observations in each interval):
    • interval specific event-indicators ${\delta }_{ij}$, where ${\delta }_{ij}=1$ if subject $i$ experienced an event in interval $j$ (i.e. $t_i\in ({\kappa }_{j-1},{\kappa }_j]$ and $T_i<C_i$) and ${\delta }_{ij}=0$ else
    • offsets $o_{ij}=\log \left(t_{ij}\right)$, where $t_{ij}=min(t_i-{\kappa }_{j-1},{\kappa }_j-{\kappa }_{j-1})$ is the time subject $i$ spent in interval $j$
  • in the $j^{th}$ interval $({\kappa }_{j-1},{\kappa }_j]$ estimate a piece-wise constant hazard rate $\lambda \left(t\right)={\lambda }_j\forall t\in ({\kappa }_{j-1},{\kappa }_j]$ (more intervals lead to better approximation)

• See Holford 1980, Laird 1981, Friedman 1982, Whitehead 1982

Application (Results)

Application (Results)

Example:

• $z=(5\times C2,6\times C3)$
• $z^{C2}=(1,1,1,1,1,0,0,0,0,0,0)$, $z^{C3}=(0,0,0,0,0,1,1,1,1,1,1)$
• $g(z,{\tilde{t}}_j=18.5)=g(z^{C2},18.5)+g(z^{C3},18.5)\approx -0.57$

$\to$ Risk reduction of $\exp (-0.57)=0.57$ compared to subject with $11\times C1$ nutrition (c.p)

Application (Results)

These bivariate surfaces are difficult to interpret as

• they must be interpreted with respect to a subject who received $C1$ nutrition on all 11 days of nutrition protocol

• partial effects $h_{C2}(t,t_e)$ and $h_{C3}(t,t_e)$ can both contribute to the cumulative effect, depending on the specific nutrition profile

• for these reasons, we prefer to analyze and interpret estimated hazard ratios between hypothetical patients with different clinically relevant exposure histories ( $z_1$ and $z_2$)

$$e_j=\frac{\lambda \left({\tilde{t}}_j|z_2\right)}{\lambda \left({\tilde{t}}_j|z_1\right)}$$

Application (Results)

We compare the following nutrition profiles:

Application (Results)

$\to$ Complete, mildly hypocaloric nutrition reduces risk of mortality compared to a complete, severely hypocaloric nutrition (Comparison B)

$\to$ No further risk reduction when moving from mildly hypocaloric to partial or complete near target nutrition (Comparisons E, F)

$\to$ Sensitivity analyses (Imputation of missing protocols, lag/lead specification, penalty structure, ...) show no substantive deviation from main results

Limitations (and outlook)

• currently, $t_{\text{lag}}$ and $t_{\text{lead}}$ must be specified a priori $\to$ would be nice if the lag-lead window could be selected data-driven (e.g. Obermeier et al., 2015)

• we assume that patients released from hospital survived until the end of the follow-up ( $t=30$ ). Sensitivity analysis with hospital discharge as censoring event do not change the results $\to$ Competing risks model for outcomes hospital discharge and death would be preferable

• Modeling and interpretation of TDCs always difficult, especially if exogeneity is unclear, e.g.

  • although nutrition is provided by hospital staff, amount provided might depend on patients' health status
  • more recent values provide better confounder adjustment but may also be fully indicative of the outcome (indication bias)

• Model choice becomes difficult when all effects are potentially non-linear and/or non-linearly time-varying (boosting ad double-penalty procedures promising)

Links and Acknowledgments

• Talk is based on two publications
  • Andreas Bender, Andreas Groll, and Fabian Scheipl. 2018. “A Generalized Additive Model Approach to Time-to-Event Analysis.” Statistical Modelling. https://doi.org/10.1177/1471082X17748083.
  • Andreas Bender, Fabian Scheipl, Wolfgang Hartl, Andrew G Day, Helmut Küchenhoff; "Penalized estimation of complex, non-linear exposure-lag-response associations", Biostatistics, , kxy003, https://doi.org/10.1093/biostatistics/kxy003

• pammtools: Package for Piece-wise exponential Additive Mixed Models (in development)

• Slides created via Yihui Xie's R package xaringan with (modified) Metropolis theme

• All graphics have been created using Hadley Whickham's ggplot2

• Models are estimated using Simon Wood's mgcv

• Web: adibender.netlify.com

• Social:

References

• Friedman, Michael. “Piecewise Exponential Models for Survival Data with Covariates.” The Annals of Statistics 10, no. 1 (1982): 101–113.
• Gasparrini, Antonio. “Modeling Exposure–lag–response Associations with Distributed Lag Non-Linear Models.” Statistics in Medicine 33, no. 5 (February 28, 2014): 881–99. https://doi.org/10.1002/sim.5963.
• Gasparrini, Antonio, Fabian Scheipl, Ben Armstrong, and Michael G. Kenward. “A Penalized Framework for Distributed Lag Non-Linear Models.” Biometrics, January 1, 2017. https://doi.org/10.1111/biom.12645.
• Holford, Theodore R. “The Analysis of Rates and of Survivorship Using Log-Linear Models.” Biometrics 36, no. 2 (1980): 299–305. https://doi.org/10.2307/2529982.
• Laird, Nan, and Donald Olivier. “Covariance Analysis of Censored Survival Data Using Log-Linear Analysis Techniques.” Journal of the American Statistical Association 76, no. 374 (1981): 231–240. https://doi.org/10.2307/2287816.
• Marra, Giampiero, and Simon N. Wood. “Coverage Properties of Confidence Intervals for Generalized Additive Model Components.” Scandinavian Journal of Statistics 39, no. 1 (March 1, 2012): 53–74. https://doi.org/10.1111/j.1467-9469.2011.00760.x.
• Sylvestre, Marie-Pierre, and Michal Abrahamowicz. “Flexible Modeling of the Cumulative Effects of Time-Dependent Exposures on the Hazard.” Statistics in Medicine 28, no. 27 (2009): 3437–3453. https://doi.org/10.1002/sim.3701.

References

• Whitehead, John. “Fitting Cox’s Regression Model to Survival Data Using GLIM.” Journal of the Royal Statistical Society. Series C (Applied Statistics) 29, no. 3 (1980): 268–75. https://doi.org/10.2307/2346901.
• Wood, Simon N. Generalized Additive Models: An Introduction with R. Boca Raton and FL: Chapman & Hall/CRC, 2006.
• Wood, Simon N. “Low-Rank Scale-Invariant Tensor Product Smooths for Generalized Additive Mixed Models.” Biometrics 62, no. 4 (December 1, 2006): 1025–36. https://doi.org/10.1111/j.1541-0420.2006.00574.x.
• Wood, Simon N. “Fast Stable Restricted Maximum Likelihood and Marginal Likelihood Estimation of Semiparametric Generalized Linear Models.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73, no. 1 (2011): 3–36. https://doi.org/10.1111/j.1467-9868.2010.00749.x.
• Wood, Simon N. “On P-Values for Smooth Components of an Extended Generalized Additive Model.” Biometrika 100, no. 1 (March 1, 2013): 221–28. https://doi.org/10.1093/biomet/ass048.
• Wood, Simon N., Fabian Scheipl, and Julian J. Faraway. “Straightforward Intermediate Rank Tensor Product Smoothing in Mixed Models.” Statistics and Computing, 2012. https://doi.org/10.1007/s11222-012-9314-z.

References

• Wickham, Hadley. Ggplot2: Elegant Graphics for Data Analysis. 2nd ed. 2016. New York, NY: Springer, 2016.
• Yihui Xie (2017). xaringan: Presentation Ninja. R package version 0.4.4. https://github.com/yihui/xaringan
• Hadley Wickham, Romain Francois, Lionel Henry and Kirill Müller (2017). dplyr: A Grammar of Data Manipulation. R package version 0.7.4. https://CRAN.R-project.org/package=dplyr

Caloric Adequacy

• caloric intake = calories from EN + PN + PF

• caloric adequacy (CA):
  $CA\left(\%\right)=\text{caloric intake}/\text{prescribed calories}\cdot 100$

• discretized caloric adequacy (in 3 categories):

  • $C1$: $0\%\le CA<30\%$ and no OI
  • $C2$:
    • $30\%\le CA<70\%$ and no OI or
    • $0\%\le CA<30\%$ and additional OI
  • $C3$:
    • $CA\ge 70\%$ or
    • $30\%\le CA<70\%$ and additional OI