Survival data with repeated systolic blood pressure measurements

Error types Likelihood Response Covariate with error Other covariate(s)
Classical, missing values Weibull survival survival time sbp1, sbp2 smoke, age, diabetes, sex
library(inlamemi)
library(INLA)
library(ggplot2)

This example shows how to fit a Weibull survival model to describe the influence systolic blood pressure (SBP) has on survival. The model is the same as in Skarstein et al. (2023), but just using inlamemi rather than plain INLA.

We assume there to be some measurement error in the SBP measurements. For some of the patients we have repeated measurements, but not for all, and for some patients both of the measurements are even missing. Therefore we are dealing with both classical measurement error and missing data in this case.

For the main model of interest, we have the formula \[ \eta_i = \beta_0 + \beta_{\texttt{sbp}} \texttt{sbp}_i + \beta_{\texttt{sex}} \texttt{sex}_i + \beta_{\texttt{age}} \texttt{age}_i + \beta_{\texttt{smoke}} \texttt{smoke}_i + \beta_{\texttt{diabetes}} \texttt{diabetes}_i. \] The error models for the repeated SBP measurement are \[ \begin{align} \texttt{sbp}^1_i = \texttt{sbp}_i + u_i^{1}, \\ \texttt{sbp}^2_i = \texttt{sbp}_i + u_i^{2}, \end{align} \] and the imputation model for sbp is \[ \texttt{sbp}_i = \alpha_0 + \alpha_{\texttt{sex}} \texttt{sex}_i + \alpha_{\texttt{age}} \texttt{age}_i + \alpha_{\texttt{smoke}} \texttt{smoke}_i + \alpha_{\texttt{diabetes}} \texttt{diabetes}_i. \] We begin by specifying the necessary priors:

# Priors for measurement error variance and true x-value
prior.prec.u <- c(0.5, 0.5) # Gamma(0.5, 0.5) (same as Keogh&Bartlett)
prior.prec.x <- c(0.5, 0.5) # Gamma(0.5, 0.5) (same as K&B)
prec.u <- 2.8
prec.x <- 1

# Prior for shape parameter of the Weibull survival model
prior.exp <- 0.01 # Gamma(1, 0.001) ~ Exp(0.001) (INLA sets prior on theta, r~Exp(0.1*theta))
exp.init <- 1.4

And then we fit the model itself. Let me point out some of the things that are special for this model:

survival_model <- fit_inlamemi(
  formula_moi = inla.surv(t, d) ~ sbp + age + smoke + sex + diabetes,
  formula_imp = sbp ~ age + smoke + sex + diabetes,
  family_moi = "weibull.surv",
  data = nhanes_survival,
  error_type = c("classical", "missing"),
  repeated_observations = TRUE,
  control.family = list(
    # Prior for main model of interest (moi)
    list(hyper = list(alpha = list(param = prior.exp,
                                   initial = log(exp.init),
                                   fixed = FALSE))),
    # Prior for error model
    list(hyper = list(prec = list(initial = log(prec.u),
                                  param = prior.prec.u,
                                  fixed = FALSE))),
    # Prior for imputation model
    list(hyper = list(prec = list(initial = log(prec.x),
                                  param = prior.prec.x,
                                  fixed = FALSE)))),
  prior.beta.error = c(0, 1/1000), # Prior for beta.sbp
  control.predictor=list(link=3)) # To specify that for the missing values, we use the third link function ("gaussian", from the imputation model) to predict them.
summary(survival_model)
#> Formula for model of interest: 
#> inla.surv(t, d) ~ sbp + age + smoke + sex + diabetes
#> 
#> Formula for imputation model: 
#> sbp ~ age + smoke + sex + diabetes
#> 
#> Error types: 
#> [1] "classical" "missing"  
#> 
#> Fixed effects for model of interest: 
#>                     mean         sd 0.025quant   0.5quant 0.975quant       mode
#> beta.0        -5.4540938 0.13791214 -5.7178638 -5.4543016 -5.1898494 -5.4543046
#> beta.age       0.9054675 0.04986843  0.8076100  0.9055121  1.0030731  0.9055131
#> beta.smoke     0.2700320 0.08453202  0.1042709  0.2700324  0.4357911  0.2700324
#> beta.sex       0.4454314 0.07880388  0.2909171  0.4454266  0.5999728  0.4454266
#> beta.diabetes  0.5919710 0.09293525  0.4096883  0.5919869  0.7741637  0.5919871
#> 
#> Coefficient for variable with measurement error and/or missingness: 
#>               mean         sd 0.025quant  0.5quant 0.975quant      mode
#> beta.sbp 0.1132246 0.04863355 0.01739928 0.1132526  0.2088873 0.1133684
#> 
#> Fixed effects for imputation model: 
#>                            mean         sd  0.025quant     0.5quant 0.975quant         mode
#> alpha.sbp.0         0.006414738 0.04311842 -0.07814269  0.006414675 0.09097253  0.006414675
#> alpha.sbp.age       0.322557254 0.02956392  0.26458304  0.322556469 0.38053594  0.322556467
#> alpha.sbp.smoke     0.004837677 0.05062415 -0.09443964  0.004837878 0.10411385  0.004837878
#> alpha.sbp.sex      -0.061813667 0.04701748 -0.15401659 -0.061814003 0.03039117 -0.061814004
#> alpha.sbp.diabetes  0.137938583 0.06224094  0.01588378  0.137937484 0.25999966  0.137937482
#> 
#> Model hyperparameters (apart from beta.sbp): 
#>                                       mean         sd 0.025quant 0.5quant 0.975quant     mode
#> Precision for sbp classical model 1.376960 0.04080456  1.2978097 1.376549   1.458454 1.376149
#> Precision for sbp imp model       2.838733 0.25689460  2.3577082 2.830223   3.368168 2.819663
#> alpha parameter for weibullsurv   1.027331 0.04959684  0.9343127 1.025741   1.129552 1.021715
plot(survival_model, plot_intercepts = FALSE, plot_imp = FALSE)
plot of chunk unnamed-chunk-5
plot of chunk unnamed-chunk-5