WAIC and LOO-CV for inference including submodels

[Tohid Jafari Koshki]

Hi all,

I want to use loo package to obtain WAIC for model comparison.

I am running a joint survival-longitudinal model written in Stan (https://groups.google.com/forum/#!category-topic/stan-users/general/1j2D34irLus). I used 3 likelihoods, one for longitudinal part or submodel and 2 for censored and uncensored data (as it is done this way in stan).

for longitudinal part: y~normal, E(y)=XB+Ub,   b is random effect and U is the coeff;

for survival: T~weibull, lambda=XB+JU, J is joining parameter.

I extract logliks for each part and then go to use loo for each of them and obtain 3 WAICs.

My question is how can I calculate single WAIC for whole model as I cannot use a single likelihood for whole dataset?

Thanks.

[Jonah Gabry]

I haven't fit any joint survival longitudinal models myself and am not very familiar with the intricacies of these models. Can you explain why the log likelihood is not just the sum of the three likelihood components you mention? If by submodels you mean you have a model for the outcome and then models also for some or all parameters in the model for the outcome, then you just have a hierarchical model and you can use the log likelihood calculated at the observation level of the model. But that's not necessarily true if that's not what submodel means here. I think Aki has a lot of experience with similar survival models so maybe he can give some more helpful advice on this particular point.

Something else relevant for loo is what is considered to be the unit left out. Is it a single measurement? A single person (so possible multiple measurements)? Etc.

Jonah

[Aki Vehtari]

Following the link you gave, I found the model but not the generated quantities block for computing likelihood terms. Please post that, too, as it would make it easier to answer your question.

Btw. we recommend using PSIS-LOO instead of WAIC, as it has better diagnostics and smaller bias.

Aki

[Tohid Jafari Koshki]

Thanks for comments. I am not sure if it is possible to put likelihoods in a single equation. I explain the model in more detail. Complete models and discussion can be found at "Xu GUO and Bradley P. CARLIN, Separate and Joint Modeling of Longitudinal and Event Time Data Using Standard Computer Packages, The American Statistician, Vol. 58, No. 1" and BUGS code in http://www.biostat.umn.edu/~brad/software.html.

For longitudinal part we use a normal distribution for each observed measurement in time

y(it)~normal(mu(i), s)

mu is modeled as mixed effects model: mu=XB+Ub  b is random effect and U is the coeff;

Then the relationship between variable y and time to event is linked using sharing components b as:

T~weibull, lambda=XB+JU, J is joining parameter.

Then in Stan (and BUGS) we write two models one for longitudinal data and another for survival.
This gives 2 DIC in BUGS each for one part.

1- I think in the article they have summed up both DICs. Could we sum 2 DICs from submodel to obtain total DIC or am I wrong?

2- How can I compute single WAIC (or similar) quantity in stan for model above?

This question is present for a simple survival model in stan where two likelihoods must be written for censored and uncensored data (as done in mice example).

I generate logliks using:

  generated quantities{

 vector [nuncens] linpreduncens;

      vector [ncens] linpredcens;

      vector [ncens] etacens;

      vector [nuncens] etauncens;

      vector [ntot] logliklong;

     vector [nuncens] loglikuncens;

     vector [ncens] loglikcens;

//generate linear predictors for longitudinal part

for (i in 1:ntot)  

           linpredlong[i]=xlong[i]*betalong+uraneff[1,subj[i]]+uraneff[2,subj[i]]*t[i];

// linpred for uncensored obs

for (i in 1:nuncens){                  

           linpreduncens[i]=xuncens[i]*betasurv;

           etauncens[i]=linpreduncens[i]+join0*uraneff[1,subj[i]]+join1*uraneff[2,subj[i]];}

// linpred for censored obs

for (i in 1:ncens){

      linpredcens[i]=xcens[i]*betasurv;

      etacens[i]=linpredcens[i]+join0*uraneff[1,subj[i]]+join1*uraneff[2,subj[i]];}

for (i in 1:ntot)

       logliklong[i]=normal_lpdf(y[i]|linpredlong[i],sigmaeps);

for (i in 1:nuncens)

           loglikuncens[i]=weibull_lpdf(tuncens[i]|alpha, exp(-etauncens[i] / alpha));

for (i in 1:ncens)

   loglikcens[i]=weibull_lpdf(tcens2[i]|alpha, exp(-etacens[i] / alpha)/tcens[i]);

  }

Many thanks and sorry if long .

[Aki Vehtari]

You can combine these likelihoods. I assume nuncens+ncens=ntot and for each y you have a time for event or censoring? Then you would probably like to sum log likelihoods for the corresponding observations and have combined log likelihood for ntot pair of observations. You have also other options but this would be most logical for me.

Aki

[Tohid Jafari Koshki]

Thank you for this. Actually I have nsubj as the total number of individuals and ntot is the total observations over nsubj as each person has multiple observation and each subject or individual fails or is censored.

Then you recommend I combine logliks of a typical person as LOGLIK=loglik(cens or event)+Sum(logliks of longitudinal observations for that person). Have I understood correctly? Then go further and use PSIS-LOO to calculate model fit.

Is there any reference for this to be cited? and/or for the other options you referred? 

Thanks for the valuable help.

[andre.p...@googlemail.com]

Tohid,

the loo is additive to the number of "observations" and requires the log-likelihood for each observation and run (iter, chains) in a matrix form.

From STAN point of view you could calculate the LL in generated quantities, but you could do it also somewhere else.

After this executing loo (waic) will give you the looic (waic) the, the lower the better. If you have different number of

observations, you could also calculate the average of looic, this comes with a danger however.

There is also the series of looic in the object and the calculation of the SE. I figured out, the LL is mostly not normal distributed.

In the documentation somewhere they note, that the number of observations mostly are big enough. Well...

Andre

[Tohid Jafari Koshki]

Thanks Andre,

Yes the number of observations differs for each subject. I wonder if obtaining LL for all observations on all subjects (ntot) and putting this matrix to loo is relevant to obtain fit statistics (waic or psis-loo) of longitudinal part? Then add waic of survival part to this to obtain overall fit waic. Is this method sensible or not? And/or is this method valid for DIC?

Previous studies have used CPO for checking model fit. But none of them has an argument on this issue.

Thank you and others who involved in this thread.

[andre.p...@googlemail.com]

Tohid,

LOO (WAIC), BIC, AIC, ... can be considered as weighted log-likelihood.

You can sum up LOOIC. The best method in my experience (although time consuming)

is leave on out-cross validation with a final validation test set.

As final warning is never good to just consider a number and say model A is better

than B, because of a single weighted LL value.

The loo method in rstan provides also a structure for each observation (sample) object.

I do know nothing about DIC. 

Andre

[Aki Vehtari]

Both following options are valid: a) leave-one-observation-out, b) leave-one-subject-out. Which one to choose depends on which prediction task is important for you a) predicting new observations for the same subjects or b) predicting new observations for new subjects. In a) you don't sum any logliks and in b) you sum all logliks of one subject (it doesn't matter if there are different number of observations and loglik terms for each subject). The case b) is however more difficult and if there are many observations for a subject it's more likely that PSIS-LOO fails (waic is even more likely to fail). Luckily loo package will report khat values so that you can check whether case b works too.

CPO is just a different name for loo log predictive density given by loo package.

DIC is relevant only if you are interested in plug-in predictive performance using the posterior mean point estimate for parameters instead predictive performance obtained by integrating over the parameters. Since you are using Stan, I assume you are integrating over the parameters and thus should not use DIC.

Aki