Survival data with repeated systolic blood pressure measurements
Source:vignettes/nhanes_survival.Rmd
nhanes_survival.Rmd| Error types | Likelihood | Response | Covariate with error | Other covariate(s) |
|---|---|---|---|---|
| Classical, missing values | Weibull survival | survival time |
sbp1, sbp2
|
smoke, age,
diabetes, sex
|
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.4And then we fit the model itself. Let me point out some of the things that are special for this model:
-
inla.surv(): Since we have a survival model, the response of the model isinla.surv(t, d). In this case,tis the survival time, anddis the censoring indicator, indicating whether the patient was still alive at the end of the study period, or whether the patient had actually passed away. -
control.family: Another thing to note is that since thefit_inlamemifunction does not have arguments for passing the prior for the shape parameter of the Weibull survival model toinla, we instead need to write out the wholecontrol.familyargument and pass this tofit_inlamemi. If you are not used to R-INLA this may look a bit strange, but this is simply how the priors for the three different levels of the model are passed toinla. So as you can see it is a list of three lists, and each of these layers corresponds to one model layer, so the first one is for the main model of interest, the second one is for the error model, and the third layer is for the imputation model. -
repeated_observations: Since we have repeated measurements for SBP, we need to set this argument toTRUE.
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)