Problems with LOO-CV in Survival

Isa_Nazar · September 25, 2021, 2:14pm

Dear all,
I want to compare the following two models using the loo package in the RSTAN. This post is a summary of my previous post (LOO-CV for Survival model - #3 by Isa_Nazar).

The following results are obtained from fitting two models under the same conditions (iter=4000, warmup=2000, thin=10, N=34744, n_grid=1001, num=12 ,province=31). It should be noted that with increasing iterations, there was no change in the results, including the value of LOOIC. My purpose is here the model assessment and choosing the best fitting model on my real data.

My survival model:

1. My survival model871×549 158 KB

Computed from 200 by 34744 log-likelihood matrix

  Estimate    SE
elpd_loo -116668.3 216.2
p_loo 63.9 0.8
looic 233358.6 432.4
Monte Carlo SE of elpd_loo is 0.6.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 34728 100.0% 157
(0.5, 0.7] (ok)

Competing model: Spatial survival model with BYM (CAR prior + Random effect)

2. Spatial survival with BYM871×549 148 KB

Computed from 200 by 34744 log-likelihood matrix

  Estimate    SE
elpd_loo -12552.8 147.8
p_loo 38.8 0.7
looic 25105.5 295.6
Monte Carlo SE of elpd_loo is 0.4.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 34738 100.0% 180
(0.5, 0.7] (ok) 6 0.0% 195
(0.7, 1] (bad) 0 0.0%
(1, Inf) (very bad) 0 0.0%
All Pareto k estimates are ok (k < 0.7).
See help(‘pareto-k-diagnostic’) for details.

The only difference between the two models is log(a[county[i]]), which has negative values, but there is a big difference between the loo values of the two models, and I was really confused. My question is whether I calculated the loo values correctly or not? Can anyone help me with this problem?
Any suggestions would be greatly appreciated!
Thank you

My stan codes used to fit the two spatial survival models are as follows:

1:My survival model:

data {
  int <lower=0> province;
  int <lower=0> n_grid;
  int <lower=0> N;
  int <lower=0> num;
  int <lower=0> number;
  int <lower =0 , upper = province> county[N];
  vector <lower=0> [province] p;
  row_vector[n_grid] kernel[province];
  matrix [N,number] respon;
}
parameters{
  vector [num] beta;
  vector <lower=0> [province] sigma1;
  real <lower=0> sigma;
  vector [n_grid] x;
  vector <lower=0> [province] r_county; 
     }
transformed parameters{
  vector [num] time_ratio;
  vector [N] lambdaa;
  real <lower=0> alpha = 1 / sigma;
  vector [n_grid] exp_x;
  vector[province] a; // `a` values between 0 and 1
  time_ratio  = exp(beta);
  exp_x = exp(x);
  lambdaa = exp (((-1 * respon[,7:18] * beta) / sigma)+log(r_county [county]));
  // }
     { // `z` not part of output
  vector[province] landa;
  for (k in 1 : province) {
        landa[k] = kernel[k] * exp_x * p[k];
     }
  a = landa / sum(landa); // assign `a` between 0 and 1
     }
     }
model{
   vector [N] s_1;
   vector [N] s_2;
   vector [N] s_new1;
   vector [N] s_new2;
   target += normal_lpdf(x|0,5);
   target += normal_lpdf(beta | 0, 1000);
   target += gamma_lpdf(alpha | 0.001, 0.001);
   target += normal_lpdf(r_county| 0, sigma1);
   target += student_t_lpdf(sigma1| 3,0,1);
for (i in 1:N) {
   s_1[i]= 1- ((lambdaa [i]*(respon[i,2]+0.01)^alpha)/ (1+(lambdaa [i]*(respon[i,2]+0.01)^alpha)));
   s_2[i]= 1- ((lambdaa [i]*(respon[i,3]+0.01)^alpha)/ (1+(lambdaa [i]*(respon[i,3]+0.01)^alpha)));
   s_new1[i] = a[county[i]] * s_1[i];
   s_new2[i] = a[county[i]] * s_2[i];
   target += log(s_1[i] - s_2[i]) + log(a[county[i]]);
   }
   }
generated quantities{
   vector [N] log_lik;
   vector [N] s_1;
   vector [N] s_2;
   vector [N] s_new1;
   vector [N] s_new2;
for (i in 1:N){
   s_1[i]= 1- ((lambdaa [i]*(respon[i,2]+0.01)^alpha)/ (1+(lambdaa [i]*(respon[i,2]+0.01)^alpha)));
   s_2[i]= 1- ((lambdaa [i]*(respon[i,3]+0.01)^alpha)/ (1+(lambdaa [i]*(respon[i,3]+0.01)^alpha)));
   s_new1[i] = a[county[i]] * s_1[i];
   s_new2[i] = a[county[i]] * s_2[i];
   log_lik [i]= log(s_1[i] - s_2[i]) + log(a[county[i]]);
  }
  }

Competing model: Spatial survival model with BYM (CAR prior + Random effect)

data {
  int <lower=0> province;
  int <lower=0> N;
  int <lower=0> num;
  int <lower=0> number;
  int <lower =0 , upper = province> county[N];
  matrix [N,number] respon;
  matrix[province,province] H;  // area adjacency matrix row standardized (ea row sums to 1)
  matrix[province,province] Delta;  // diagonal matrix in BYM 1991, sfstd model B1
  matrix[province,province] Ones; // diag matrix of ones for identity
  real phi; // phi = 1/maxeigval(Delta*H*Delta)
     }
transformed data {
  vector[province] zeros;  // zeros for mean of MVN specification
for (i in 1:province)
  zeros[i] = 0.0;
     }
parameters{
  vector [num] beta;
  real <lower=0> sigma;
  vector[province] w1;  // random effect (CAR)
  real<lower=1e-5> prec_sre;  // precision of CAR
  vector[province] theta;       // heterogeneous effects
     }
transformed parameters{
  vector [num] time_ratio;
  vector [N] lambdaa;
  real <lower=0> alpha = 1 / sigma;
  real <lower=0> tau_sq;
  time_ratio  = exp(beta);
  tau_sq = 1.0 / (prec_sre * prec_sre);
  lambdaa = exp ((-respon[,7:18] * beta - w1[county] - theta[county]) / sigma);
     }
model{
   vector [N] s_1;
   vector [N] s_2;
   matrix[province,province] Sigma_sre;  // covariance matrix for CAR
   target += normal_lpdf(beta | 0, 1000);
   target +=gamma_lpdf(alpha |0.001,0.001);
   prec_sre ~ gamma(0.001,0.001);    // precision of spatial re
// spatial re
   Sigma_sre = (tau_sq * inverse(Ones - phi * H)) * Delta;
   w1 ~ multi_normal(zeros, Sigma_sre);
   target += normal_lpdf(theta| 0, 1);
for (i in 1:N) {
   s_1[i]= (1- ((lambdaa [i]*(respon[i,2]+0.01)^alpha)/ (1+(lambdaa [i]*(respon[i,2]+0.01)^alpha))));
   s_2[i]= (1- ((lambdaa [i]*(respon[i,3]+0.01)^alpha)/ (1+(lambdaa [i]*(respon[i,3]+0.01)^alpha))));
   target += log(s_1[i] - s_2[i]);
    }
    }
generated quantities{
   vector [N] log_lik;
   vector [N] s_1;
   vector [N] s_2;
  for (i in 1:N){
    s_1[i]= (1- ((lambdaa [i]*(respon[i,2]+0.01)^alpha)/ (1+(lambdaa [i]*(respon[i,2]+0.01)^alpha))));
    s_2[i]= (1- ((lambdaa [i]*(respon[i,3]+0.01)^alpha)/ (1+(lambdaa [i]*(respon[i,3]+0.01)^alpha))));
    log_lik [i]= log(s_1[i] - s_2[i]);
    }
    }
.```

martinmodrak · September 30, 2021, 12:55pm

Hi,
the models are a bit too complex for me to understand quickly, so I’ll presume there are no bugs (you should definitely check, e.g. by testing if you can recover parameters from simulated data). I’ll also note that my understanding of LOO is relatively surface-level, but since nobody else answered, I will give it a try.

One way to try to understand why LOO values are very different is to note that LOO works with the posterior predictive distribution. Large differences in LOO-CV should thus usually correspond to different posterior predictions for new data. Do the predicitons of the models (for some new unseen data point) differ substantially? If not, but you get large differences in LOO-CV, I would suspect a bug somewhere.

I really don’t understand what log(a[county[i]]) is doing in the likelihood, but maybe the problem really lies there - maybe the models differ in what parameters you are conditioning on when computing the log likelihood? But that’s just a very wild guess…

Best of luck with your model!

Isa_Nazar · October 2, 2021, 3:41pm

Thanks a ton for your reply and help. Here, the only difference between the likelihood of two models is log(a[county[i]]) that it is added in the log-likelihood function to account for the spatial dependency.
Is it possible to get different results by comparing models using two different GOF criteria, LOO & DIC? I mean, a model is selected in terms of LOO values, but its DIC value is more.

Thank you in advance!

jsocolar · October 2, 2021, 3:53pm

Are you sure this is achieving what you intend? In your transformed parameters block you note in a comment that a is constrained to take values between 0 and 1, which means that its logarithm is guaranteed to be negative. If this is truly the only difference between the likelihoods, then it doesn’t surprise me that LOO would dramatically prefer the model that doesn’t have a giant negative term tacked onto the likelihood.

jsocolar · October 2, 2021, 4:02pm

It’s certainly possible (otherwise we wouldn’t need all the awesome work that’s gone into LOO and we could all just use DIC!). But if it’s true that the only difference in the likelihoods here is a big negative term, then I don’t think there’s any way that DIC will give you different results here.

Isa_Nazar · October 2, 2021, 4:03pm

Thank you so much for your quick reply and assistance. Yes, these are probability values that have been estimated and their logarithm is negative and the only difference between the two models is the same. Do you have any suggestions to compare the two models in RSTAN?

Thank you

Isa_Nazar · October 2, 2021, 4:26pm

Apart from the log(a[county[i]]) values, the covariates were linked differently to the scale parameter in the two models as below:

1:My survival model:
lambdaa = exp (((-1 * respon[,7:18] * beta) / sigma)+log(r_county [county]));

Competing model: Spatial survival model with BYM (CAR prior + Random effect)
lambdaa = exp ((-respon[,7:18] * beta - w1[county] - theta[county]) / sigma);

Here, w1 is a spatial random effect with Car prior, and theta is an ordinary random effect.

Thank you

avehtari · October 6, 2021, 6:09pm

By a quick look I didn’t see anything wrong. As p_loo for both models is quite small compared to the number of observations, you don’t need LOO to see the difference between the models. It would be good to look at the actual predictions and some expected loss in scale that is easier for you to interpret.

Isa_Nazar · October 7, 2021, 1:27pm

Thanks a ton, Prof. Vehtari for your reply and clarification of my confusion. I think I need to change my pointwise log-likelihood in the competing model as follows because it’s just the likelihood of the data (conditional on the parameters, beta, sigma,sigma1, and r_county). any advice and suggestions will be greatly appreciated.
‘’’
generated quantities{
vector [N] log_lik;
vector [N] s_1;
vector [N] s_2;
// vector [N] s_new1;
//vector [N] s_new2;
for (i in 1:N){
s_1[i]= 1- ((lambdaa [i] (respon[i,2]+0.01)^alpha)/ (1+(lambdaa [i] (respon[i,2]+0.01)^alpha)));
s_2[i]= 1- ((lambdaa [i] (respon[i,3]+0.01)^alpha)/ (1+(lambdaa [i] (respon[i,3]+0.01)^alpha)));
//s_new1[i] = a[county[i]] * s_1[i];
// s_new2[i] = a[county[i]] * s_2[i];
log_lik [i]= log(s_1[i] - s_2[i]) ;
}
}
‘’’.

Thank you in advance!

Warm Regards,
Eisa

Topic		Replies	Views
LOO-CV for Survival model Modeling rstan , fitting-issues , loo	3	556	September 23, 2021
I used stan_jm with different option for assoc like etavalue, etaauc etc. How can I compare these models, any fitness information I can get from these model? rstanarm	5	843	January 22, 2021
LOOIC diagnostics with Stanfit object Modeling loo	5	638	October 30, 2020
Errors when calculating PSIS-LOO CV in loo version 2.2.0 Modeling fitting-issues , loo	7	1236	January 2, 2020
Model checking & comparison using loo vs loo_compare Modeling	8	1843	January 4, 2024

Problems with LOO-CV in Survival

Related topics