How can I get quantiles of predictions in INLA

[Samar Deen]

Following the FAQ Help Topic: "How can I get predictions using INLA", I was able to predict posterior means. How can I get the 0.025quant, and 0.975quant to make the plot given below?

I would eventually want to do out of sample predictions. But that will be a next step.

Data Files attached:

data4<-read.csv("Falldata4.csv")

data4$y1<-exp(data4$y)-1

Fluke.adj <- paste(getwd(),"/Fluke.graph",sep="")
g = inla.read.graph("Fluke.graph")

formula.zip<- y1 ~ BotTemp+BotSalin+
  f(ID.area,model="bym", graph=g) +
  f(ID.year,model="rw2") +
  f(ID.year1,model="iid")+
  f(ID.area.year,model="iid")

link = rep(NA, nrow(data4))
link[which(is.na(data4$y1))] = 1
result.pred = inla(formula.zip,
                   data = data4, offset = log(data4$E+1),
                   family = "zeroinflatedpoisson1",
                   control.predictor=list(compute=T))

n = nrow(data4)
fitted.values.mean = numeric(n)
for(i in 1:n) {
  if (is.na(data4$y1[i])) {
    if (FALSE) {
      ## either like this, which is slower...
      marg = inla.marginal.transform(
        function(x) exp(x)/(1+exp(x)),
        result.pred$marginals.fitted.values[[i]] )
      fitted.values.mean[i] = inla.expectation(
        function(x) x, marg)
    } else {
      ## or like this,  which is much faster...
      fitted.values.mean[i] = inla.emarginal(
        function(x) exp(x)/(1 +exp(x)),
        result.pred$marginals.fitted.values[[i]])
    }
  } else {
    fitted.values.mean[i] = result.pred$summary.fitted.values[i,"mean"]
  }
}

fitted.vals<-as.data.frame(fitted.values.mean[1:n])

#Plot Results

data.pred<-cbind(data4,fitted.vals)

data.pred$Year<-data4$ID.year

plot(supsmu(data.pred$Year, data.pred$fitted.values.mean[1:n]), type = "l", col = "red", ylim = c(-0.1,2), ylab = "Fluke (ALB Units)", xlab = "Years(1991-2020)")
lines(supsmu(data.pred$Year, data.pred$`0.025quant`), col = "red", lwd=0.5, lty=2)
lines(supsmu(data.pred$Year, data.pred$`0.975quant`), col = "red", lwd=0.5, lty=2)
lines(supsmu(data.pred$Year, data.pred$y1), col = "black", lwd=0.5)
abline(v=24)

[Haakon Bakka]

If I understand you correctly, the difficulty is with integrating out the likelihood distribution / randomness.

Perhaps this will help you: https://haakonbakka.bitbucket.io/btopic112.html#adding-observation-likelihood-randomness

- Here I do integration via random samples (Monte Carlo integration I think it's called)

PS: This example also discuss out of sample prediction.

Kind regards,

Haakon Bakka

--
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-group+unsub...@googlegroups.com.
To post to this group, send email to r-inla-discussion-group@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.

[Samar Deen]

Hi Haakon,

Thank you for your email. I have trying diffferent variations of the code you've provided to get the quantiles for out of sample predictions, but I havent had much luck. This is what I've done so far. Any pointers or help is most appreciated.

Basically I want to show how the outcome variable (Number of Fish) will vary, once I change the Bottom Temperature at a rate of 0.23 celcius per year, I want to demonstrate how the prediction interval also changes in out of sample predictions.

#This merges the samples to the main dataset - Can I use these anyway?

n.samples = 1000

samples = inla.posterior.sample(n.samples, result = result.pred)

df<-as.data.frame(samples[[1]]$latent[1:2592]) #one sample

samples.pred = lapply(samples, function(x) x$latent[1:2592]) #get all the samples

with.samples<-cbind(samples.pred,data4)

#I am not sure how I can get out of sample predictions with this?

our.experiment = list(BotTemp = 6, BotSalin = 36, ID.area = 7, ID.year =24,
                      ID.year1 = 24, ID.area.year = 2491 )
predictor.our.experiment<-NULL

for (nr in 1:length(samples)){

s = samples[[nr]]$latent
f.bt = s[5450, , drop=F] * our.experiment$BotTemp
f.bs = s[5451, , drop=F] * our.experiment$BotSalin
f.id.area = s[2593, , drop=F]
f.id.year = s[2809, , drop=F]
f.id.year1 = s[2857, , drop=F]
f.id.area.year = s[5449, , drop=F]
int = s[5450, , drop=F]

predictor.our.experiment[nr]  = drop(f.bt + f.bs + f.id.area + f.id.year+ f.id.year1+f.id.area.year+int)
}

samples.link = inla.link.invlogit(predictor.our.experiment)
plot(density(samples.link))

To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion-group+unsubscr...@googlegroups.com.

To post to this group, send email to r-inla-discussion-group@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.

--

Samar Deen

Summer Research Fellow,  Atkinson Center for a Sustainable Future
PhD Candidate, Department of Natural Resources, Cornell University

M.P.A. Cornell University, 2011

Phone: +1 607 793 4200
Skype ID: samarnasiruddeen

[Haakon Bakka]

Hi Samar,

I am not entirely sure how to reply to this. My best reply would be to re-state the things I wrote in the webpage.

To help yourself construct the right samples, you want to write down the mathematical expression for the linear predictor.

Your code looks roughly reasonable. Each experiment will need its own for-loop, where you have a different Bottom Temperature for each "experiment".

I guess I can do it for you if you want, but then you must send me all the code and the data so I can run it.

PS: There is a different way to solve this problem that may be easier in some cases (and impossible in other cases): Namely using the marginals.fitted.values, transform them and add likelihood noise.

Kind regards,

Haakon

[Haakon Bakka]

One small thing though:

"#This merges the samples to the main dataset - Can I use these anyway?"

I don't see why you try to do this. The posterior samples and the data are two very different things.