Prior predictive checks in survival modeling with a mixture prior

Modeling
, , ,
1

Hi all,

I recently posted for help to build a survival model using a mixture prior. I was able to run this model without issue, code below:

data {
  int<lower=0> K;          // num covariates

  int<lower=0> N;          // num uncensored obs
  vector[N] t;             // event time (non-strict decreasing)
  matrix[N, K] x;          // covariates for uncensored obs

  int N_c;                 // num censored obs
  vector[N_c] t_c;         // censoring time
  matrix[N_c, K] x_c;      // covariates for censored obs
}

parameters {
  vector[K] beta;          // slopes (no intercept)
}

model {
  // mixture of 0.4*N(0,4) + 0.6*N(0,0.24)
  for (k in 1:K) {
    target += log_mix(0.4, normal_lpdf(beta[k] | 0, 4),
                           normal_lpdf(beta[k] | 0, 0.24));
  };

  vector[N] log_theta = x * beta;
  vector[N_c] log_theta_c = x_c * beta;

  real log_denom = log_sum_exp(log_theta_c);

  for (n in 1:N) {
    log_denom = log_sum_exp(log_denom, log_theta[n]);
    target += log_theta[n] - log_denom;   // log likelihood
  }

}

I was also hoping to generate draws from this mixture for prior predictive checks (and to make sure the distribution is behaving as expected).

A minimal working example in R (using dnorm rather than rnorm) looks like this:

# the x axis of the plot
x <- seq(from = -15, to = 15, length.out = 2000)

# the y axes of the plot
## Y ~ N(0, 4)
y_n04 <- dnorm(x, 0, 4)

## Y ~ N(0, 0.24)
y_n0024 <- dnorm(x, 0, 0.24)

## Y ~ 0.4*N(0, 4) + 0.6*N(0, 0.24)
y_mix <- y_n04*0.4 + y_n0024*0.6

# plot the distributions
plot(x, y_n04, type = "l", ylim = c(0, 2))
lines(x, y_n0024)
lines(x, y_mix)

I was unable to replicate this type of generated prior check within the generated quantities{} block in Stan. I tried several different combinations of normal_rng(), normal_lpdf(), and log_mix() without success (either getting errors or PRNG values that did not replicate the expected mixture distribution.

Any direction would be appreciated and thanks again, all, for your time and patience!

Best,
Mazin

2

This question must have been an issue with my beginner status in Stan. Commenting out the rest of the model outputs the expected prior distribution when plotted:

data {
  int<lower=0> K;          // num covariates

  int<lower=0> N;          // num uncensored obs
  vector[N] t;             // event time (non-strict decreasing)
  matrix[N, K] x;          // covariates for uncensored obs

  int N_c;                 // num censored obs
  vector[N_c] t_c;         // censoring time
  matrix[N_c, K] x_c;      // covariates for censored obs
}

parameters {
  vector[K] beta;          // slopes (no intercept)
}

model {
  // mixture of 0.4*N(0,4) + 0.6*N(0,0.24)
  for (k in 1:K) {
    target += log_mix(0.4, normal_lpdf(beta[k] | 0, 4),
                           normal_lpdf(beta[k] | 0, 0.24));
  };

  // vector[N] log_theta = x * beta;
  // vector[N_c] log_theta_c = x_c * beta;

  // real log_denom = log_sum_exp(log_theta_c);

  // for (n in 1:N) {
  //   log_denom = log_sum_exp(log_denom, log_theta[n]);
  //   target += log_theta[n] - log_denom;   // log likelihood
  // }

}

Plotting the beta output overlays well with the theoretical mixture distribution.

Thanks all!
Mazin