% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/full.R
\name{conformal.pred}
\alias{conformal.pred}
\title{Conformal prediction intervals.}
\usage{
conformal.pred(
  x,
  y,
  x0,
  train.fun,
  predict.fun,
  alpha = 0.1,
  w = NULL,
  mad.train.fun = NULL,
  mad.predict.fun = NULL,
  num.grid.pts = 100,
  grid.factor = 1.25,
  verbose = FALSE
)
}
\arguments{
\item{x}{Matrix of features, of dimension (say) n x p.}

\item{y}{Vector of responses, of length (say) n.}

\item{x0}{Matrix of features, each row being a point at which we want to
form a prediction interval, of dimension (say) n0 x p.}

\item{train.fun}{A function to perform model training, i.e., to produce an
estimator of E(Y|X), the conditional expectation of the response variable
Y given features X. Its input arguments should be x: matrix of features,
y: vector of responses, and out: the output produced by a previous call
to train.fun, at the \emph{same} features x. The function train.fun may
(optionally) leverage this returned output for efficiency purposes. See
details below.}

\item{predict.fun}{A function to perform prediction for the (mean of the)
responses at new feature values. Its input arguments should be out: output
produced by train.fun, and newx: feature values at which we want to make
predictions.}

\item{alpha}{Miscoverage level for the prediction intervals, i.e., intervals
with coverage 1-alpha are formed. Default for alpha is 0.1.}

\item{w}{Weights, in the case of covariate shift. This should be a vector of
length n+n0, giving the weights (i.e., ratio of test to training feature
densities), at each of n+n0 the training and test points. Default is NULL,
which means that we take all weights to be 1.}

\item{mad.train.fun}{A function to perform training on the absolute residuals
i.e., to produce an estimator of E(R|X) where R is the absolute residual
R = |Y - m(X)|, and m denotes the estimator produced by train.fun.
This is used to scale the conformal score, to produce a prediction interval
with varying local width. The input arguments to mad.train.fun should be
x: matrix of features, y: vector of absolute residuals, and out: the output
produced by a previous call to mad.train.fun, at the \emph{same} features
x. The function mad.train.fun may (optionally) leverage this returned
output for efficiency purposes. See details below. The default for 
mad.train.fun is NULL, which means that no training is done on the absolute
residuals, and the usual (unscaled) conformal score is used. Note that if
mad.train.fun is non-NULL, then so must be mad.predict.fun (next).}

\item{mad.predict.fun}{A function to perform prediction for the (mean of the)
absolute residuals at new feature values. Its input arguments should be
out: output produced by mad.train.fun, and newx: feature values at which we
want to make predictions. The default for mad.predict.fun is NULL, which
means that no local scaling is done for the conformal score, i.e., the
usual (unscaled) conformal score is used.}

\item{num.grid.pts}{Number of grid points used when forming the conformal
intervals (each grid point is a trial value for the interval). Default is
100.}

\item{grid.factor}{Expansion factor used to define the grid for the conformal
intervals, i.e., the grid points are taken to be equally spaced in between
-grid.factor*max(abs(y)) and grid.factor*max(abs(y)). Default is 1.25. In
this case (and with exchangeable data, thus unity weights) the restriction
of the trial values to this range costs at most 1/(n+1) in coverage. See
details below.}

\item{verbose}{Should intermediate progress be printed out? Default is FALSE.}
}
\value{
A list with the following components: pred, lo, up, fit. The first
  three are matrices of dimension n0 x m, and the last is a matrix of
  dimension n x m. Recall that n0 is the number of rows of x0, and m is the
  number of tuning parameter values internal to predict.fun. In a sense, each
  of the m columns really corresponds to a different prediction function;
  see details below. Hence, the rows of the matrices pred, lo, up give
  the predicted value, and lower and upper confidence limits (from conformal
  inference), respectively, for the response at the n0 points given in
  x0. The rows of fit give the fitted values for the n points given in x.
}
\description{
Compute prediction intervals using conformal inference.
}
\details{
For concreteness, suppose that we want to use the predictions from
  forward stepwise regression at steps 1 through 5 in the path. In this case,
  there are m = 5 internal tuning parameter values to predict.fun, in the
  notation used above, and each of the returned matrices pred, lo, and up will
  have 5 rows (one corresponding to each step of the forward stepwise path).
  The code is structured in this way so that we may defined a single pair of
  functions train.fun and predict.fun, over a set of m = 5 tuning parameter
  values, instead of calling the conformal function separately m = 5 times.

  The third arugment to train.fun, as explained above, is the output produced
  by a previous call to train.fun, but importantly, at the \emph{same}
  features x. The function train.fun may (optionally) leverage this returned
  output for efficiency purposes. Here is how it is used, in this function: 
  when successively querying several trial values for the prediction
  interval, denoted, say, y0[j], j = 1,2,3,..., we set the out argument in
  train.fun to be the result of calling train.fun at the previous trial value
  y0[j-1]. Note of course that the function train.fun can also choose to 
  ignore the returned output out, and most default training functions made
  available in this package will do so, for simplicity. An exception is the
  training function produced by \code{\link{lm.funs}}. This will use this 
  output efficiently: the first time it trains by regressing onto a matrix of
  features x, it computes an appropriate Cholesky factorization; each
  successive time it is asked to train by regressing onto the same matrix x,
  it simply uses this Cholesky factorization (rather than recomputing one).
  The analogous explanation and discussion here applies to the out argument
  used by mad.train.fun.

  If the data (training and test) are assumed to be exchangeable, the basic
  assumption underlying conformal prediction, then the probability that a new
  response value will lie outside of [-max(abs(y)), max(abs(y))], where y is
  the vector of training responses, is 1/(n+1).  Thus the restriction of the
  trials values to [-grid.factor*max(abs(y)), grid.factor*max(abs(y))], for
  all choices grid.factor >= 1, will lead to a loss in coverage of at most
  1/(n+1). This was also noted in "Trimmed Conformal Prediction for
  High-Dimensional Models" by Chen, Wang, Ha, Barber (2016) (who use this
  basic fact as motivation for proposing more refined trimming methods).
}
\examples{
## Lasso: use a fixed sequence of 100 lambda values

# Generate some example training data
set.seed(33)
n = 100; p = 120; s = 10
x = matrix(rnorm(n*p),n,p)
beta = c(rnorm(s),rep(0,p-s)) 
y = x \%*\% beta + rnorm(n)

# Generate some example test data
n0 = 50
x0 = matrix(rnorm(n0*p),n0,p)
y0 = x0 \%*\% beta + rnorm(n0)

# Grab a fixed lambda sequence from one call to glmnet, then
# define lasso training and prediction functions
if (!require("glmnet",quietly=TRUE)) {
  stop("Package glmnet not installed (required for this example)!")
}
out.gnet = glmnet(x,y,nlambda=100,lambda.min.ratio=1e-3)
lambda = out.gnet$lambda
funs = lasso.funs(lambda=lambda)

# Conformal inference, and jacknife and split conformal versions
out.conf = conformal.pred(x, y, x0, alpha=0.1,
  train.fun=funs$train, predict.fun=funs$predict, verb=TRUE)

out.jack = conformal.pred.jack(x, y, x0, alpha=0.1,
  train.fun=funs$train, predict.fun=funs$predict)

out.split = conformal.pred.split(x, y, x0, alpha=0.1, seed=0,
  train.fun=funs$train, predict.fun=funs$predict)

y0.mat = matrix(rep(y0,ncol(out.conf$lo)),nrow=n0)
cov.conf = colMeans(out.conf$lo <= y0.mat & y0.mat <= out.conf$up)
len.conf = colMeans(out.conf$up - out.conf$lo)
err.conf = colMeans((y0.mat - out.conf$pred)^2)

cov.jack = colMeans(out.jack$lo <= y0.mat & y0.mat <= out.jack$up)
len.jack = colMeans(out.jack$up - out.jack$lo)
err.jack = colMeans((y0.mat - out.jack$pred)^2)

cov.split = colMeans(out.split$lo <= y0.mat & y0.mat <= out.split$up)
len.split = colMeans(out.split$up - out.split$lo)
err.split = colMeans((y0.mat - out.split$pred)^2)

# Compare to parametric intervals from oracle linear regression
lm.orac = lm(y~x[,1:s])
out.orac = predict(lm.orac,newdata=list(x=x0[,1:s]),
  interval="predict", level=0.9)

cov.orac = mean(out.orac[,"lwr"] <= y0 & y0 <= out.orac[,"upr"])
len.orac = mean(out.orac[,"upr"] - out.orac[,"lwr"])
err.orac = mean((y0 - out.orac[,"fit"])^2)
  
# Plot average coverage 
plot(log(lambda),cov.conf,type="o",pch=20,ylim=c(0,1),
     xlab="log(lambda)",ylab="Avg coverage",
     main=paste0("Conformal + lasso (fixed lambda sequence):",
       "\nAverage coverage"))
points(log(lambda),cov.jack,type="o",pch=20,col=3)
points(log(lambda),cov.split,type="o",pch=20,col=4)
abline(h=cov.orac,lty=2,col=2)
legend("bottomleft",col=c(1,3,4,2),lty=c(1,1,1,2),
       legend=c("Conformal","Jackknife conformal",
         "Split conformal","Oracle"))

# Plot average length
plot(log(lambda),len.conf,type="o",pch=20,
     ylim=range(len.conf,len.jack,len.split,len.orac),
     xlab="log(lambda)",ylab="Avg length",
     main=paste0("Conformal + lasso (fixed lambda sequence):",
       "\nAverage length"))
points(log(lambda),len.jack,type="o",pch=20,col=3)
points(log(lambda),len.split,type="o",pch=20,col=4)
abline(h=len.orac,lty=2,col=2)
legend("topleft",col=c(1,3,4,2),lty=c(1,1,1,2),
       legend=c("Conformal","Jackknife conformal",
         "Split conformal","Oracle"))

# Plot test error
plot(log(lambda),err.conf,type="o",pch=20,
     ylim=range(err.conf,err.jack,err.split,err.orac),
     xlab="log(lambda)",ylab="Test error",
     main=paste0("Conformal + lasso (fixed lambda sequence):",
       "\nTest error"))
points(log(lambda),err.jack,type="o",pch=20,col=3)
points(log(lambda),err.split,type="o",pch=20,col=4)
abline(h=err.orac,lty=2,col=2)
legend("topleft",col=c(1,3,4,2),lty=c(1,1,1,2),
       legend=c("Conformal","Jackknife conformal",
         "Split conformal","Oracle"))

}
\references{
See "Algorithmic Learning in a Random World" by Vovk, Gammerman,
  Shafer (2005) as the definitive reference for conformal prediction; see
  also "Distribution-Free Predictive Inference for Regression" by Lei,
  G'Sell, Rinaldo, Tibshirani, Wasserman (2018) for another description; and
  "Conformal Prediction Under Covariate Shift" by Barber, Candes, Ramdas,
  Tibshirani (2019) for the weighted extension.
}
\seealso{
\code{\link{conformal.pred.jack}},
  \code{\link{conformal.pred.split}}, \code{\link{conformal.pred.roo}}
}
