Introduction to conformalbayes

Source:vignettes/conformalbayes.Rmd

conformalbayes.Rmd

Theconformalbayes package provides functions toconstruct finite-sample calibrated predictive intervals for Bayesianmodels, following the approach inBarber et al. (2021). Thebasic idea is a natural one: use cross-validated residuals to estimatehow large predictive intervals need to be, on average.

Suppose we have a heavy-tailed dataset.

library(conformalbayes)library(rstanarm)library(ggplot2)sim_data=function(n=50){x=rnorm(n)y=3-2*x+rt(n, df=2)data.frame(x=x, y=y)}d_fit=sim_data()ggplot(d_fit,aes(x,y))+geom_point()+geom_smooth(method=lm, formula=y~x)

We can fit a linear regression to the data, but it won’t give usaccurate uncertainty quantification in our predictions.

# fit the modelm=stan_glm(y~x, data=d_fit, chains=1, refresh=0)d_test=sim_data(2000)interv_model=predictive_interval(m, newdata=d_test, prob=0.50)# are the points coveredcovered_model=with(d_test,interv_model[,1]<=y&y<=interv_model[,2])ggplot(d_test,aes(x,y, color=covered_model, group=1))+geom_point(size=0.4)+geom_linerange(aes(ymin=interv_model[,1],                       ymax=interv_model[,2]), alpha=0.4)+labs(color="Covered?")+geom_smooth(method=lm, formula=y~x, color="black")

In fact, the 50% intervals over-cover, with a coverage rate of 69.8%,since the fat tails of the error terms pulls the estimate of theresidual standard deviation too high.

While a posterior predictive check could uncover this discrepancy,leading us to fit a more flexible model, we can take another approachinstead. The first step is to callloo_conformal(), whichcomputes leave-one-out cross-validation weights and residuals for use ingenerating more accurate predictive intervals.

m=loo_conformal(m)print(m)#> stan_glm#>  family:       gaussian [identity]#>  formula:      y ~ x#>  observations: 50#>  predictors:   2#> ------#>             Median MAD_SD#> (Intercept)  2.9    0.3#> x           -1.9    0.3#>#> Auxiliary parameter(s):#>       Median MAD_SD#> sigma 2.1    0.2#>#> ------#> * For help interpreting the printed output see ?print.stanreg#> * For info on the priors used see ?prior_summary.stanreg#> (conformalbayes enabled, with estimated CI inflation factor 0.81 and#> transformation pair function (x) x and function (x) x)

Theloo_conformal() returns the same fitted model, justwith a thin wrapping layer that contains the leave-one-outcross-validation information. You can see at the bottom of the outputthatconformalbayes estimates that correctly-sizedpredictive intervals are only 81% of the size of the model-basedpredictive intervals.

To actually generate predictive intervals, we usepredictive_interval(), just like normal:

interv_jack=predictive_interval(m, newdata=d_test, prob=0.50)# are the points coveredcovered_jack=with(d_test,interv_jack[,1]<=y&y<=interv_jack[,2])ggplot(d_test,aes(x,y, color=covered_jack, group=1))+geom_point(size=0.4)+geom_linerange(aes(ymin=interv_jack[,1],                       ymax=interv_jack[,2]), alpha=0.4)+labs(color="Covered?")+geom_smooth(method=lm, formula=y~x, color="black")

Indeed, the coverage rate for these jackknife conformal intervals is49.2%, as we would expect.

The conformal version ofpredictive_interval() doescontain two extra options:plus andlocal.Whenplus=TRUE, the function will generatejackknife+ intervals, which have a theoretical coverageguarantee. These can be computationally intensive, so by default theyare only generated when the number of fit and prediction data points isless than 500. In practice, non-plus jackknife intervalsgenerally perform just as well asjackknife+ intervals.Whenlocal=TRUE (the default), the function will generateintervals whose widths are proportional to the underlying model-basedpredictive intervals. So if your model accounts for heteroskedasticity,or produces narrow intervals in areas of covariate space with manyobservations (like a linear model),local=TRUE will producemore sensible intervals. The overall conformal performance guaranteesare unaffected.