DanWeinberger
commented
3 years ago Cleaner code:
set.seed=(123) n = 100 alpha=3 beta =2 x= rnorm(n, 0,0.2)
y= rpois(n, exp(alpha+x*beta) )
r = inla(y ~ 1 + x, data = data.frame(y, z), control.compute = list(config=TRUE), family = "nbinomial") r.samples = inla.posterior.sample(101, r)
pred.interval.func <- function(sample.ds, dist=c('nb','poisson')){ mu1 <- exp(sample.ds$latent[grep('Predictor', row.names(sample.ds$latent))]) if(dist=='nb'){ nb.size1 = sample.ds$hyperpar pred <- replicate(102, rnbinom(n=length(mu1), mu=mu1, size=nb.size1), simplify = 'array') }else{ pred <- replicate(102, rpois(n=length(mu1), lambda=mu1), simplify = 'array') } return(pred) }
samps <- sapply(r.samples,pred.interval.func,dist='poisson', simplify='array') samps.q <- t(apply(samps, 1, quantile, probs=c(0.025,0.5, 0.975)))
par(mfrow=c(1,1))
plot(x, y) arrows(x0=x, y0=samps.q[,1],y1=samps.q[,3], length=0)
Intervals given by INLA seem to not account for observation likelihood. Need to add in additional sampling step to get prediction interval: set.seed=(123) n = 100 alpha=1 beta =2 x= rnorm(n, 0,0.2)
y= rpois(n, exp(alpha+x*beta) )
r = inla(y ~ 1 + x, data = data.frame(y, z), control.compute = list(config=TRUE), family = "poisson") r.samples = inla.posterior.sample(10^3, r)
post.labels<-dimnames(r.samples[[1]]$latent)[[1]]
posterior.samples<- sapply(r.samples, '[[', 'latent') #should be 'joint' instead of 'latent' https://www.ucl.ac.uk/population-health-sciences/sites/population-health-sciences/files/inla_baio.pdf
preds.select<-grep('Predictor',post.labels )
posterior.preds <- exp(((posterior.samples[preds.select,]))) #lambda
pred.q <- t(apply(posterior.preds,1,quantile, probs=c(0.025,0.5,0.975)))
pois.pred.func <- function(){ matrix(rpois(length(posterior.preds), lambda=posterior.preds), nrow=nrow(posterior.preds)) } posterior_pred_pois <- pbreplicate(101,pois.pred.func(), simplify = 'array' )
posterior.preds.pois.q <- t(apply(posterior.preds.pois,1,quantile, probs=c(0.025,0.5,0.975)))
cbind(posterior.preds.pois.q,y)
par(mfrow=c(1,1)) plot(posterior.preds.pois.q[,2], y) abline(a=0, b=1) arrows(x0=posterior.preds.pois.q[,1],, x1=posterior.preds.pois.q[,3],y0=y, length=0)