CPO and PIT for predictions

[Noah Silverman]

Apologies for re-posting.  My last post was flooded by spammers, and I don't think anyone here read it.

I'm fitting a model, and setting some values to NA, to test out-of-sample predictive performance.

Normally, I use PIT and CPO to help diagnose model fit.  However, for NA values INLA doesn't provide those (That's expected, as there's no ground-truth to use for calculation)

Is there a way to reconcile this?  Perhaps retroactively provide the true values, and then get the CPO and PIT?

Thank You

[Helpdesk]

On Mon, 2022-09-05 at 02:13 -0700, Noah Silverman wrote:
> I'm fitting a model, and setting some values to NA, to test out-of-
> sample predictive performance.
>
> Normally, I use PIT and CPO to help diagnose model fit.  However, for
> NA values INLA doesn't provide those (That's expected, as there's no
> ground-truth to use for calculation)
>

in this case, the posterior for that linear predictor is what you need
in order to compute the prediction of any response. but with NA there is
no way to do that, unless you construct a fake response yourself.

> Is there a way to reconcile this?  Perhaps retroactively provide the
> true values, and then get the CPO and PIT?
>
> Thank You

> --
> You received this message because you are subscribed to the Google
> Groups "R-inla discussion group" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to r-inla-discussion...@googlegroups.com.
> To view this discussion on the web, visit
> https://groups.google.com/d/msgid/r-inla-discussion-group/a04fe578-9f26-489b-8da4-e37a3493821fn%40googlegroups.com
> .

--
Håvard Rue
he...@r-inla.org

[Noah Silverman]

That makes sense.

But, We DO HAVE the ground truth values for this data.  We wanted to test the model performance on unseen data.  So, we simulated the unseen, but changing some values to NA, then re-running INLA to see how accurately the model can capture those data.

Is there a function we can extract from INLA, or some manual way to calculate CPO and PIT for specific data points?
We have:
- Inla predictions (distributions) for the NA values
- True values (in separate data structure)

We can easily calculate RMSE, but since INLA predictions are actually a distributions, something that captures that, like CPO or PIT would seem a better metric.

Any ideas?

Thanks!

[Helpdesk]

it should be the same to add the true values for the NA's as the cpo
computes

\pi(y_i | y_{-i})

so that this data-point, or true values, will be removed when the CPO
for data_i is computed.

[Kyle Manning]

In the case that the data has a temporal aspect, are CPO/ PIT still alright metrics to use? 

Does a leave-one-out strategy still make sense for such temporally dependent cases? 

Another question - if data is naturally clustered, is it ok to use CPO/ PIT/ leave one out metrics?

For example, if we measure test results of students across different countries, and group students by class:

if we include random effects on a class level, can we still use such metrics?

If I imagine fitting a model leaving only one student from a class out, but I have information from the other students, won't results be biased

Or, does the estimation method account for such possibilities?

[Subramanian Swaminathan]

Hi team,

I am having similar problem: computing CPO and PIT for out-sample predictions. Can we use the " inla.posterior.sample" and "inla.posterior.sample.eval" commands to predict the response and the corresponding CPO and PIT values for the out of samples given the fitted "spatiotemporal model" , covariates (X1, X2,X3) and the actual values for the response in place of "NA" used in the fitted model ("result"). If so, in the "fun" command how the model and covariates corresponding to the out-of sample prediction to be incorporated. I was trying to use the code given in the " R-INLA help" but was not successful. I would appreciate if some one in the group provide the code with an example.

Thanks in advance

ssubra

[Helpdesk (Haavard Rue)]

it is often easier to add these new prediction points into the model
with an NA response, and then you can work it out from there directly

[Subramanian Swaminathan]

Thanks for your prompt response. I tried your suggestion incorporating NAs for the new prediction points and refitting the model. The cpo and pit values were found to be 'NA' for all these new prediction points. It appears that I need to give the actual response values for calculating the cpo and pit for these new prediction points. Refitting the model with actual response means that I am fitting the model including the new prediction points. kindly advise. 

Best regards

ssubra

[Helpdesk (Haavard Rue)]

you have to do this last step manually. for a new point 'j', then you'll
get

\pi( \eta_j | y)

so then you need to integrate this into the likelihood, like

\int \pi(y_j | \eta_j ) \pi(\eta_j | y) d\eta_j

like by approximating it with a sum

> https://groups.google.com/d/msgid/r-inla-discussion-group/34cd16c1-914c-452c-a593-7db010e3a66en%40googlegroups.com
> .

[Finn Lindgren]

But this can only be done in that way for distributions that only depend on eta, and doesn’t have any hyperparameters (theta). So for Poisson, this works, but not for Gaussian observations where the variance/precision parameter is estimated.

A solution to this is to use posterior sampling to get joint samples of eta and theta, and to use that to estimate the response level probabilities/densities for the new test observations. We’ve done in the past by computing the average of samples of e.g.

dnorm(y, mean= eta_expression, sd=(The_precision_parameter_name)^-1)

Where eta_expression and the precision parameter are compute or extracted from the posterior sample outputs. Similar can be done for any observation model, including zero-inflated models, etc, but it does require some user-side coding.

Finn

> On 14 Sep 2022, at 15:18, Helpdesk (Haavard Rue) <he...@r-inla.org> wrote:
>
> 

> To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/8e4d887a1e7eb98babd496f50df696d0282cad35.camel%40r-inla.org.

[Helpdesk (Haavard Rue)]

You're right.

If the point where you want to predict is in the model then 'Predictor'
and 'APredictor' (when using stack/spde-models) are predefined for each
sample.

[Kyle Manning]

Hi, 

I have a similar question about predictive distributions/ sampling.

First off, I have a more general question: in a setting where I have a Gaussian/ T distribution, 

my understanding is I can access the estimated precision parameter using:

samples = inla.posterior.sample(n, model, num.threads= <>)

Then this gives me the estimated precision:

samples[[i]]$hyperpar[1]

I guess I was wishing to have different posterior predictive distribution to have different variances, but it seems that I will get the same precision for every observation.

It would be desirable to have different distributions/ variances, is there a way to implement this? Or, am I misunderstanding and are the precisions not the same across observations?

Finn, you mention that you use posterior sampling. I am still not quite sure how to go about getting the posterior predictive distribution.

Is it correct to say that functions like smarginal() dmarginal(), etc. are not reflective of the true predictive distribution, but rather the distribution of the linear predictor? 

I have read on this site and Julian Faraways book about using an extra random effect to estimate these distributions, is that the preferred method?

Or, it seems like there is an alternative method - what do you mean by using the average of samples of dnorm() ? 

That was a lot of questions - let me know your thoughts!

Thanks,

Kyle

[Håvard Rue]

if you're interested in sampling hyperparameters, then use

inla.hyperpar.sample

as 'inla.posterior.sample' only use the integration points (which is why
it seems like the hyperparameters are the same, which they are, by
design...)

> https://groups.google.com/d/msgid/r-inla-discussion-group/2aa1af66-ad81-4fb8-b47c-1160bdfd7181n%40googlegroups.com
> .

--
Håvard Rue
hr...@r-inla.org

[Subramanian Swaminathan]

Hi Finn and Rue,

Many thanks for all your valuable suggestions!

I have fitted 9 additive models (for areal data) using icar, or bym2 or lcar priors for space and iid or RW1 or RW2 for time.  I have been trying to calculate the PIT values for assessing the predictive performance of the model. As suggested by r-inla help  team, I tried to calculate PIT from the fitted models. I  have drawn the PIT histograms for 3 models (icar, bym2 & lcar priors) with iid prior for time. All the 3 PIT histograms seems to be same.  I do not know if my PIT calculations are correct? Though I have used the 'theta'  for 'size' of the 'negbin' in the calculation, I do not know how to plug-in the other 'theta' values corresponding to space and time in the calculation. I would appreciate your kind advice on what changes are required in the attached code. 

Further, I am getting warnings as  "There were 50 or more warnings (use warnings() to see the first 50)". Whether the warnings indicate that the predictions are not available for some of the records.

The code, PIT histograms and the log files are attached for your kind information. 

Best regards

ssubra

[Helpdesk (Haavard Rue)]

yes, the plots looks a but 'suspicious'...

are you computing the PIT values like done with

control.compute=list(cpo=TRUE)

?

in this case, it PIT values are already computed. (You may want to
correct with discrete data, since its <= Y, hence you can subtract 1/2
the cpo value

> https://groups.google.com/d/msgid/r-inla-discussion-group/b7f0b654-cf7f-4d95-9011-7ed087fe0f52n%40googlegroups.com
> .

[Subramanian Swaminathan]

HI Prof Rue,

Thanks for your response!

While fitting the models I have used "control.compute=list(cpo=TRUE)'.  I have fitted the models to training data, keeping 'NA' for the response variable in the test data.  Now in order to have an assessment of the predictive performance of the model, I am trying to calcuate the PITs for the test data, for which there are no PIT values available in the output, as I have used 'NA' while model fitting. Therefore, I have adapted the modified version of the codes for discrete distribution as per Czado  et al. Biometrics, 2009. In this, the cumulative probabilities and probabilities of the negative binomial model are calculated by plugging in the estimates of 'mu' and 'size' of the NB. The estimates of the covariates are obtained from the model fitted with training data. I was able to calcualte the PITs for a fixed effect model, by plugging in the estimates of the covariates in "mu" and "size'"(please see the code in PIT.R file attached in my earlier mail). However, I do not know how the 'mu' to be estimated for additive models (spatial, temporal random effects) and spatiotemporal interaction models. For your kind information, before initiating the calculations, I loaded the .RData image file to get relevant estimates for the PIT calculation.

I would be grateful if you could suggest any other methods of calculation.

Best regards

ssubra

[Helpdesk (Haavard Rue)]

I attach a simple example, hopefully this will be helpful

> https://groups.google.com/d/msgid/r-inla-discussion-group/08a1ce07-abe6-4539-ae77-fc007d62ea1dn%40googlegroups.com
> .

[Subramanian Swaminathan]

Hi Prof Rue,

Appreciate your valuable suggestions with simple example (runme4.R). That works fine. However, I do not understand whether 'p' values (if I am right these are PITs)  from the posterior samples are for the 'test data' OR for the 'entire data'. When I checked the 'p' values, I found it was available for the 'entire data set'. I am interested in assessing the predictive performance (in terms of PIT histograms) of different models corresponding to 'test data'. Kindly let me know, whether I can draw the PIT histogram corresponding to the 'test data' and if so how can I filter the P values for the test data. The relevant codes from your simple example are given below for your quick reference:

## compute PIT values, the easy (approximative) way. we compuare with the cpo/pit values which

## are not the same (since data are remove), but its about the same with a lot of data (and a
## simple model like here).

xx <- inla.posterior.sample(1000, r)
fun <- function(E, Y) {
    ## theta[1] is 'size',  this must be verified from 'summary(r)'
    size <- theta[1]
    prob <- size / (size + E * exp(Predictor))
    return (pnbinom(Y, size = size, prob = prob))
}
## here we insert complete data
p <- rowMeans(inla.posterior.sample.eval(fun, xx, E = E, Y = y))
plot(p, r$cpo$pit, pch = 19)
abline(a = 0, b = 1, lwd = 3, col = "blue")
hist(p,  breaks=10,main="test data",xlab="PIT")

Once gain thanks a lot for your great help!

ssubra

[Helpdesk (Haavard Rue)]

On Thu, 2022-10-06 at 03:14 -0700, Subramanian Swaminathan wrote:
> ## compute PIT values, the easy (approximative) way. we compuare with
> the cpo/pit values which
> ## are not the same (since data are remove), but its about the same
> with a lot of data (and a
> ## simple model like here).
>
> xx <- inla.posterior.sample(1000, r)
> fun <- function(E, Y) {
>     ## theta[1] is 'size',  this must be verified from 'summary(r)'
>     size <- theta[1]
>     prob <- size / (size + E * exp(Predictor))
>     return (pnbinom(Y, size = size, prob = prob))
> }
> ## here we insert complete data
> p <- rowMeans(inla.posterior.sample.eval(fun, xx, E = E, Y = y))
> plot(p, r$cpo$pit, pch = 19)
> abline(a = 0, b = 1, lwd = 3, col = "blue")
> hist(p,  breaks=10,main="test data",xlab="PIT")

its easiest to compute it like this, for the complete data, argument 'Y
= y', but you know which to predict, so you can just extract that part
of 'p'

--
Håvard Rue
he...@r-inla.org

[Subramanian Swaminathan]

Hi Prof Rue,

As suggested, I have filtered the PITs for the testing data and prepared the PIT histograms. It works fine. I really appreciate your timely help always! GREAT!

Best regards

ssubra

[Subramanian Swaminathan]

Dear Prof. Rue,

In continuation of your suggestions, I have extracted the PIT values corresponding to unobserved data. Similarly, I would like to calculate the CPO values for the unobserved data.  Could you please guide me with some simple example as you had given for PIT? or would it be possible to get both CPO and PIT in one go?

 Appreciate your kind advice always,

ssubra

[Helpdesk (Haavard Rue)]

The CPO could be computed with the PIT, as PIT is

PIT = Prob(Y.new <= y.obs)

while

CPO = Prob(Y.new = y.obs)

[Subramanian Swaminathan]

Dear Prof. Rue,

I got it, thanks a lot. 

One clarification, can you use 'dnbinom' (negbin likelihood) function to calculate the 'CPO' for the test data. I have cross checked the calculation for the training data with that of CPO obtained from the model output and it matches. I would be grateful if you could clarify this too. 

Once again thank you for your kind response.

Regards

ssubra

[Helpdesk (Haavard Rue)]

if you want the CPO then you use dnbinom, yes.

[Subramanian Swaminathan]

Thanks a lot! Appreciate your valuable timely response. 

Happy that i am improving in R -INLA with your kind immense helps in time. 

Best wishes

ssubra