REIN, July 2023

Author

damiano.cerasuolo @ unicaen.fr

Population

Event/Populations
Entire population	1 549 / 45 026
ADPKD(-)	1 505 / 40 961
ADPKD(+)	44 / 2 560

Power Generalised Weibull Distribution

The PGW distribution (Nikulin and Haghighi, 2009) is a three-parameter distribution.

The Weibull distribution is a continuous probability distribution that can fit an extensive range of distribution shapes. Like the normal distribution, the Weibull distribution is unimodal and describes probabilities associated with continuous data. However, unlike the normal distribution, it can also model skewed data. In fact, its extreme flexibility allows it to model both left- and right-skewed data.

Its survival function can be written as follows:

\[ S(t; \sigma; \nu; \gamma) = exp \{ 1- [1 +(\frac{t}{\sigma})^\nu ]^{1/\gamma} \} \]

Parameters used in the PGW are:

sigma : scale parameter
nu : shape parameter
gamma : shape parameter
t : positive argument
p : probability (0,1)
n : number of simulations

The PGW adds a second shape parameter to the Weibull distribution.

Time in the data set

Quantiles of follow-up time (days)

  0%  10%  20%  30%  40%  50%  60%  70%  80%  90% 100% 
  10  364  394  722  727  734 1087 1094 1453 1470 1825

Histogram of follow-up time

The PGW allows to model non-normal distribution (skewed unimodal) of time-to-event, but not to model interval censoring. To get around the interval censoring, a data transformation could be applied.

After the exclusion of 6 219 patients with a follow-up time < 10 days, censoring occurs as follows:

year
1	4 107
2	11 932
3	7 434
4	7 2921
5	7 525
	38 922

A normal distribution is attributed to each interval. Mean and standard deviation are arbitrary tough.

qX <- rnorm_bounded(38922, mean = Z, sd = P, min = V, max = W)

rrdbfr <- rrdbfr %>% 
  mutate(timenom = case_when(
    timecl == "R" ~ qX, 
    [...]
  ))

A rnorm_bounded version of normal distribution function is used since for each time lapse a min and max values are used. The distribution does not appear to be normal, multi-modal nor skewed. Neverthless, the PGW can be applied.

Variable matrix

# Design matrices
des <- as.matrix(cbind(scale(rrdbfr$age), rrdbfr$sex, rrdbfr$apkd01))
colnames(des) <- cbind("std age", "sex", "apkd01")
des_t <- as.matrix(cbind(scale(rrdbfr$age)))

p0 <- dim(des_t)[2]
p1 <- dim(des)[2]

# Required quantities
status <- as.logical(rrdbfr$status)
times <- as.vector(rrdbfr$timenom)/365.24 # in years
times.obs <- times[status]

Adaptive Metropolis within Gibbs sampler

chain <- adaptMetropGibbs(ltd=log.post, starting=inits, accept.rate=0.44, batch=n.batch, 
                  batch.length=batch.length, report=50000, verbose=FALSE)$p.theta.samples

Trace plots of posteriors samples (after thinning and burn-in period)

burn <- 5000
thin <- 50

Burn-in describes the practice of throwing away some iterations at the beginning of an MCMC run. The burn-in notion says that (1) you start somewhere (i.e., at x); (2) you run the Markov chain for n steps; (3) you throw away all the data (no output). This is the burn-in period. After the burn-in you run normally, using each iterate in your MCMC calculations.¹

To sum up things, burn-in (or warm-up), by removing early values of the Markov chain, eliminates the impact of the starting value.

There are some alternatives to burn-in, though not explored here.

Thinning has nothing to do with Bayesian inference, but everything to do with computer-based pseudo-random simulation. By subsampling or thinning out the Markov chain the correlation between the simulations of the posterior distibution (which is the final purpose of Markov chain generated via MCMC algorithms) is reduced - as the thinning interval grows.

However, it has nothing to do with convergence of the Markov chain to the stationary distribution since it is a post-processing of the simulated Markov chain.

It makes sense once the chain is (approximately) stationary. ²

Histograms of posteriors samples (after thinning and burn-in period)

Result : marginal survival function using NMC³ ⁴ ⁵ posterior samples

Posterior sampling

Posterior distribution is the result of the prior information and the likelihood. Sampling from this distribution allows approximating posterior expected values, quantiles, marginal densities, etc.

With NMC = 1000, the posterior survival distribution is:

plot(survfit(Surv(times, status) ~ 1, data = rrdbfr), 
     xlab = "Time", 
     ylab = "Marginal survival probability")
curve(msurv,0,14, ylim = c(0,1), lty = 1, lwd = 2, add = T, col = "blue")
legend("topright",  legend = c("GH Model", "Kaplan-Meier"), col = c("blue","black"), lwd = c(2,2) )

where msurv is calculated as follow:

msurv <- Vectorize(function(t){
  temp <- vector()
  for(j in 1:NMC){
    exp.x.alphap <- as.vector(exp(des_t%*%chainp[j,3]))
    exp.x.beta.difp <- as.vector(exp( des%*%chainp[j,(3+p0):(2+p0+p1)] - des_t%*%chainp[j,3] ))
    temp[j] <- mean(exp(- chpgw(t*exp.x.alphap, chainp[j,1], chainp[j,2], chainp[j,3])*exp.x.beta.difp)) 
  }

Questions/To do-list

Why the model is so f**ked up compared to KM estimator? Is the GH wrong? Is the KM wrong?
Some parameters clearly do not fit (par 2, i.e.). How to evaluate and improve that?
Variables to be included in the prior distribution
Model output (HR, 95%CI)

Bayes Survival: unadjusted analysis with dependent \(\gamma\) prior

With Gamma (\(\alpha\), \(\beta\)), we refer to the Gamma distribution with shape parameter \(\alpha\) and rate parameter \(\beta\).

res <- BayesSurv(df = rdb.ba, #our data frame
                 time = "time.zero", #survival/censoring times
                 event = "status.n", #status indicator
                 prior = "Dependent", #dependent Gamma prior
)

# plot bayes: hazard
PlotBayesSurv(bayes.surv.object = res,
              object = "survival")

Graphic representations

Hazard
Survival
Cumulative Hazard

Time in BayesSurvival⁶

In the analysis time = time + 0.00001. This is because survSplit does not support zero time when waiting time is 0 (in that case, start time would be < endtime).

Mixed model with `BMRS`

mod1 <- stan_glmer(status ~ apkd01 + dial + bmic + sex + age + diabetes + 
                      (1 | timeclass), data = rdb)

Model information:

random intercept: time class (see the table below)
random spline: no
prior: Gaussian

time value (days)	class
0	0
1-365	1
366-727	2
727-1143	3
1144-1509	4
≥1510	5

Model output


Model Info:
 function:     stan_glmer
 family:       gaussian [identity]
 formula:      status ~ apkd01 + dial + bmic + sex + age + diabetes + (1 | timeclass)
 algorithm:    sampling
 sample:       4000 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 36345
 groups:       timeclass (6)

Estimates:
                                           mean   sd   10%   50%   90%    mcse Rhat n_eff
(Intercept)                               1.0    0.0  1.0   1.0   1.0     0.0  1.0   778 
apkd011                                   0.0    0.0  0.0   0.0   0.0     0.0  1.0  5753
dial2                                     0.0    0.0  0.0   0.0   0.0     0.0  1.0  5901
bmic2                                     0.0    0.0  0.0   0.0   0.0     0.0  1.0  1681
bmic3                                     0.0    0.0  0.0   0.0   0.0     0.0  1.0  1665
bmic4                                     0.0    0.0  0.0   0.0   0.0     0.0  1.0  1738
sex2                                      0.0    0.0  0.0   0.0   0.0     0.0  1.0  6161
age                                       0.0    0.0  0.0   0.0   0.0     0.0  1.0  5225
diabetes1                                 0.0    0.0  0.0   0.0   0.0     0.0  1.0  5290
b[(Intercept) timeclass:0]                0.0    0.0 -0.1   0.0   0.0     0.0  1.0   859 
b[(Intercept) timeclass:1]                0.1    0.0  0.0   0.1   0.1     0.0  1.0   695
b[(Intercept) timeclass:2]                0.0    0.0  0.0   0.0   0.0     0.0  1.0   694 
b[(Intercept) timeclass:3]                0.0    0.0  0.0   0.0   0.0     0.0  1.0   693
b[(Intercept) timeclass:4]                0.0    0.0  0.0   0.0   0.0     0.0  1.0   689
b[(Intercept) timeclass:5]                0.0    0.0  0.0   0.0   0.0     0.0  1.0   712
sigma                                     0.2    0.0  0.2   0.2   0.2     0.0  1.0  3717
Sigma[timeclass:(Intercept),(Intercept)]  0.0    0.0  0.0   0.0   0.0     0.0  1.0   880 
mean_PPD                                                                  0.0  1.0  4234
log-posterior                                                             0.1  1.0   926

Fit Diagnostics:
           mean   sd   10%   50%   90%
mean_PPD   1.0    0.0  1.0   1.0   1.0

Joint model

The mean_ppd is the sample average posterior predictive distribution of the outcome variable or the distribution of possible unobserved values conditional to the observed values.

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

A simpler model

mod.a <- stan_glmer(status ~ apkd01  + (1 | timeclass), data = rdb)

It is possible to select other priors (other than Gaussian):

mod <- stan_glmer(Y ~ X1 + X2 + X3 + (1 | Z), data = data,
                  prior = student_t(df = n), prior_intercept = student_t(df = n),
                  seed = 21051986)

Image from the shinystan for the model. (launch_shinystan(mod.a).)

Joint model in `rstanarm` and `survival`⁷

Joint modelling can be broadly defined as the simultaneous estimation of two or more statistical models which traditionally would have been separately estimated.

They include:

a longitudinal mixed effects model which analyses patterns of change in an outcome variable that has been measured repeatedly over time (for example, a clinical biomarker) and
a survival or time-to-event model which analyses the time until an event of interest occurs (for example, death or disease progression).

Joint estimation of these so-called “submodels” is achieved by assuming they are correlated via individual-specific parameters (i.e. individual-level random effects).

Data sets

data type	data set name	size
longitudinal data	rl	348 972
survival data	reinONE	24 516

Model exemple

Research question

What is the research question?

stan_jm(formulaLong = age ~ cardiovasc + delyears.num + (delyears.num | ID),
                dataLong = rl,
                formulaEvent = survival::Surv(time, status) ~ 1,
                dataEvent = rlONE,
                time_var = "delyears.num",
                chains = 1, refresh = 2000, seed = 12345)

Pareto distribution error

Error in if (p$diagnostics$pareto_k > 1) { : 
  missing value where TRUE/FALSE needed

Footnotes

http://users.stat.umn.edu/~geyer/mcmc/burn.html↩︎
https://stats.stackexchange.com/questions/442714/why-does-thinning-work-in-bayesian-inference↩︎
NMC = Newtonian Monte Carlo (Arora, et al., 2020) is a second-order gradient-based Markov chain Monte Carlo (MCMC) algorithm that uses the first- and second-order gradients to propose a new value for a random variable↩︎
NMC theorizes that the posterior is shaped like a Normal distribution↩︎
https://beanmachine.org/docs/newtonian_monte_carlo/↩︎
https://cran.r-project.org/web/packages/BayesSurvival/vignettes/VignetteBayesSurvival.html↩︎
https://mc-stan.org/rstanarm/articles/jm.html↩︎

Population

Power Generalised Weibull Distribution

Time in the data set

Quantiles of follow-up time (days)

Histogram of follow-up time

Variable matrix

Adaptive Metropolis within Gibbs sampler

Trace plots of posteriors samples (after thinning and burn-in period)

Histograms of posteriors samples (after thinning and burn-in period)

Result : marginal survival function using NMC3 4 5 posterior samples

Bayes Survival: unadjusted analysis with dependent \(\gamma\) prior

Graphic representations

Mixed model with BMRS

Model output

A simpler model

Joint model in rstanarm and survival7

Data sets

Model exemple

Footnotes

Result : marginal survival function using NMC³ ⁴ ⁵ posterior samples

Mixed model with `BMRS`

Joint model in `rstanarm` and `survival`⁷