[2]Bertrand Clarke

1]Charles Schwab Financial, Street, Denver, 100190, CO, USA

2]Department of Statistics, U. Nebraska-Lincoln, 340 Hardin Hall North,
Lincoln, 68583-0963, NE, USA

# Post-model-selection prediction for GLM’s

Dean Dustin<EMAIL_ADDRESS><EMAIL_ADDRESS>[  [

###### Abstract

We give two prediction intervals (PI) for Generalized Linear Models that take
model selection uncertainty into account. The first is a straightforward
extension of asymptotic normality results and the second includes an extra
optimization that improves nominal coverage for small-to–moderate samples.
Both PI’s are wider than would be obtained without incorporating model
selection uncertyainty. We compare these two PI’s with three other PI’s. Two
are based on bootstrapping procedures and the third is based on a PI from
Bayes model averaging. We argue that for general usage either the asymptotic
normality or optimized asymptotic normality PI’s work best. In an Appendix we
extend our results to Generalized Linear Mixed Models.

###### keywords:

prediction interval, genealized linear model, post-model selection

##  1 Introduction

It is well known that linear models and their extensions - generalized linear,
linear mixed, generalized linear mixed models, and generlized mixed models -
are the workhorses of statistical analysis. Aside from formulating such
models, analysts have to chose amongst competing models usually in the same
class. Unless a model is proposed pre-experimentally, model selection from a
model list is done after the data is collected. There are numerous model
selection procedures, but regardless of which model is chosen as “best”, the
resulting model will have an associated variability inherited from the
variability in the data. How to take this variability into account properly
when making predictions is the main topic of this paper.

Common practice in many predictive contexts is to choose a model and then use
it to generate predictions. Such plug-in methods are common in many modeling
contexts. This is pragmatic but neglects taking account of the uncertainty
dure to model selection or, pehaps more commonly, variable selection. Here we
propose prediction intervals (PI’s) for generalized linear models (GLM’s) that
are modified by the model selection principle (MSP) used for variable
selection so that their nomial coverage is asymptotically correct in the limit
of sample sizes. This is important because Hong et al (2018a)) showed that
using model selection procedure procedures like Akaike’s information criterion
can result in predictive intervals with lower than nominal coverage if the
PI’s do not take the uncertainty of the MSP into account.

The post model selection inference problem has gained wide interest since the
problem was first addressed in Berk et al (2013). The so called post selection
inference (PoSI) intervals introduced in Berk et al (2013) are universally
valid for any model selection principle (MSP). However, PoSI intervals are
known to be conservative (see Leeb et al (2015)) partially because they allow
for any ad-hoc MSP to be used. The PoSI framework was used to construct
universally valid (over all MSP’s) confidence intervals for the mean of a
predictive distribution in LM’s in Bachoc et al (2019). Universally valid
confidence regions for the simultaneous inference problem are constructed
Kuchibhotla et al (2020). A different approach was proposed by Efron (2014)
that uses bootstrap intervals to address the post model selection inference
problem under a single choice of MSP.

Stine (1985) introduced bootstrapped predictive intervals in linear regression
in 1985, but these intervals did not consider uncertainty due to model
selection. Leeb (2009) introduced a model selection procedure based on cross
validation techniques and proved that using this technique, the resulting
prediction interval from the selected model is approximately valid. While this
is a seemingly strong and useful result, it holds only in the high sparsity
case with $n<<p$ as well as in the limit of large $n$. Specifically in his
Proposition 4.3 the intervals are guaranteed to be within
$1/\sqrt{n}+\epsilon$ of the nominal coverage for $0<\epsilon\leq\log(2)$.
More recently, predictive intervals based on Shorth – i.e. the shortest
interval containing a pre-specified number of values – for GLM’s and GAMs are
studied in Stine (2021). While these intervals are valid, and account for
uncertainty due to the MSP, they are not as intuitive and general as the ones
we present in Subsec. 3.

Our methodology is in contrast to the PoSI-based intervals from Berk et al
(2013) that essentially widen PI’s until the nomial coverage is achieved.
Indeed, the PoSI intervals take the pessimistic (if practical) view that model
developers will use MSP’s that are not theoretically sound. Our approach is
optimistic in that we assume a proper MSP with well known theoretical
properties will be used. This allows us to incorporate the variability from a
given generic MSP into our PI’s.

Here, we present two PI’s that account for the uncertainty of an MSP in an
intuitive manner for GLM’s. These PI’s are easy to understand and,
importantly, are easy to implement. We also present a less intuitive way to
construct a PI to give better coverage along the lines of PoSI intervals. This
PI seems work well in terms of coverage, but interpretation the interval is
difficult.

The structure of this paper is as follows. In Sec. 2 we define the notation
and setting needed for our approach. In fact, the notation incorporates much
of the intuition behind our approach. In Sec. 3, we present the main theorem
that gives a PI’s dependent on an MSP for the case of GLM’s. We then give a
finite sample improvement for use with this PI in small-to-moderate sample
settings. We also define three other intervals, two based on bootstrapping and
one from the Bayes model average. In Sec. 4 we present our simulation results.
Finally, in Sec. 5 we summarize the implications of our work. We extend our
theory to GLMM’s in the Appendices.

##  2 Notation and Setting

Throughout this paper we assume model selection and variable selection are
synonymous, and defined as followings. Let
$\mathcal{D}_{n}=\\{(y_{1},x_{1}),\ldots,(y_{n},x_{n})\\}$ where $Y_{i}=y_{i}$
is an outcome of the response variable and $x_{i}$ is a value of the
$d$-dimensional explanatory variable. We use superscripts to indicate vectors,
thus $y^{n}=(y_{1},\ldots,y_{n})^{T}$. Let $m\in\mathcal{M}$ be a candidate
model in the full collection of models $\mathcal{M}$. We define a variable
selection procedure $M=M(\mathcal{D}_{n})$ which takes the available data and
maps it to a subset of variables based on some objective function we denote
$Q$. We denote a chosen model $\hat{m}=\arg\min_{m}Q(m,\mathcal{D}_{n})$. We
think of $Q$ as an objective function such as the Akaike of Bayes information
criterion (AIC, BIC) or as a penalized loss function. For instance, for linear
models, if $Q$ is the AIC, we have

$Q_{AIC}(m,\mathcal{D}_{n})=-2\ln(p(y^{n}|X^{n}_{m},\hat{\beta}^{MLE}_{m}))+2d,$
(2.1)

where $d$ is the number of parameters that need to be estimated in model
$m\in\mathcal{M}$. This defines a function

$M:\mathcal{D}_{n}\mapsto\mathcal{M}.$

from the data $\mathcal{D}_{n}$ into the model space. That is, we think of the
variable selection procedure as a function
$M:\mathbb{R}^{n}\times\mathbb{R}^{n\times d}\mapsto\mathcal{M}$.

Other choices for $Q$ include the BIC given by

$Q_{BIC}(m,\mathcal{D}_{n})=2\ln(p(y^{n}|X^{n}_{m},\hat{\beta}^{MLE}_{m}))+d\log(n),$
(2.2)

and the general Bethel-Shumway class of information criteria defined in Bachoc
et al (1988). The BIC and the Bethel-Shumway class of information criteria are
consistent for model/variable selection.

###  2.1 Variable Selection

In the predictive context, the interpretation of a parameter in a linear model
is consistent across models: If the parameter, say, $\beta_{j}$, appears in
multiple models it always means the change in $Y$ for a unit change in $x_{j}$
holding other explanatory variables constant. This is in contrast to the
modeling view, see Berk et al (2013), that regards each parameter as an
element in a whole model. and affected by each other’s values. Thus, for us,
if $\beta_{j}=0$, this is mathematically equivalent to a model that does not
include $x_{j}$. Indeed, in practice, we often set a threshold $\eta>0$ and
and say that when $|\beta_{j}|<\eta$, we set $\beta_{j}=0$.

With this in mind, we can define $M$ as follows. Let $X=(x_{1},\ldots,x_{d})$,
and define $M=\\{\hat{\delta}_{1},\ldots,\hat{\delta}_{d}\\}$ where for
$j=1,\ldots,d$

$\hat{\delta}_{j}=\begin{cases}1\text{ if }X_{j}\text{ is selected
under}~{}$Q$\\\ 0\text{ otherwise. }\end{cases}$

For the true model we have $m_{T}=\\{\delta_{1,T},\ldots,\delta_{p,T}\\}$
where

$\delta_{j,T}=\begin{cases}1\text{ if }X_{j}\in m_{T}\\\ 0\text{ otherwise.
}\end{cases}$

Define the set $\mathcal{M}_{S}$ to be the set containing all possible
combinations $M$ can take. The cardinality of $\mathcal{M}_{S}$ is $2^{d}$ and
assuming $m_{T}$ exists, $m_{T}\in\mathcal{M}_{S}$. We write

$\beta_{T}=(\beta_{\delta_{1,T}},\ldots,\beta_{\delta_{p,T}})\text{ and
}\hat{\beta}_{M}=(\hat{\beta}_{\hat{\delta}_{1}},\ldots,\hat{\beta}_{\hat{\delta}_{p}}),$

and $\dim(\hat{\beta}_{M})=\dim(\beta_{m_{T}})$. Note that if
$\hat{\delta}_{j}=0$, then by default we set
$\hat{\beta}_{\hat{\delta_{j}}}=0$. Furthermore, we have that $\delta_{j,T}=0$
is equivalent to $\beta_{j}=0$.

We specify our target of inference as $\beta_{T}$ and write

$\beta_{m_{T}}=\left(X^{\prime}_{m_{T}}X_{m_{T}}\right)^{-1}X^{\prime}_{m_{T}}E(Y)$

in the linear models context. That is we are trying to estimate the true
parameters, regardless of the model that is chosen. This is in contrast to the
random target of inference

$\beta_{M}=\left(X^{\prime}_{M}X_{M}\right)^{-1}X^{\prime}_{M}E(Y)$

defined in Berk et al (2013). Thus, for linear models, we define the estimate
for $\beta_{j}$ as follows:

$\hat{\beta}_{\hat{M},j}=\begin{cases}\left[\left(X^{\prime}_{\hat{M}}X_{\hat{M}}\right)^{-1}X^{\prime}_{\hat{M}}y\right]_{j}\text{
if }\hat{\delta}_{j}=1\\\ 0\hskip 97.56493pt\text{ if
}\hat{\delta}_{j}=1.\end{cases}$

Now there are two steps in the process of obtaining the true model. The first
step is to estimate the $\delta_{j}$’s. In this step we want $M$ to give
$\hat{\delta}_{j}=1$ if $\delta_{j,T}=1$, however, $M$ may also give
$\hat{\delta}_{j}=1$ even if $\delta_{j,T}=0$. In this case, our definition
allows the estimate $\hat{\beta}_{\delta_{j}}$ to be zero. Thus, even if $M$
includes variables that are not in $M_{T}$ we can still estimate their
coefficients to be zero which allows $M\rightarrow m_{T}$ asymptotically (as
seen in Theorem 3.1).

###  2.2 Prediction in Generalized Linear Models

As noted, we restrict attention to GLM’s and GLMM’s. To be more precise,
suppose $Y\sim\mathcal{G}(\mu,R)$ where $\mathcal{G}$ is an exponential family
with mean $\mu$ and variance $R$. Then the pdf of $Y$ given the canonical
parameter $\theta$ is

$f(y|\theta)=e^{\frac{y\theta-b(\theta)}{a(\phi)}+c(y,\phi)}$ (2.3)

where $\phi$ is a scale parameter. In one parameter exponential families such
as Poisson or Binomial distributions, $a(\phi)=1$. From (2.3) we have the
following properties:

  * •

$E(Y|X)=\frac{\partial b(\theta)}{\partial\theta}=\mu$

  * •

$Var(Y|X)=a(\phi)\frac{\partial^{2}b(\theta)}{\partial\theta^{2}}=a(\phi)V(\mu)$

  * •

$I(\theta)=Var(\ell(\theta|y,\phi))$

Following standard GLM practice, we model the mean of $Y$ by transforming it
to a linear function of the explanatory variables. The function we use to
transform $E(Y)=\mu$ is called the link function and we denote it by
$g(\cdot)$. Note that $g(\cdot)$ is a continuous invertible function. This
gives us the linear predictor

$\eta=g(E(Y|X))=g(\mu)=X\beta$ (2.4)

and we define the inverse link function to be the inverse of $g(\cdot)$ which
is

$\mu=E(Y|X)=g^{-1}(X\beta).$ (2.5)

For nwo, we assume that $X$ is of full rank, to avoid problems with
estimability.

Note that the canonical parameter $\theta$ is a function of $\mu$ so we write
$\theta=\theta(\mu)=\theta(g^{-1}(X\beta))$. Now we can write the log-
likelihood of (2.3) as

$\ell(\beta|y,\phi)=\frac{y\left(\theta(g^{-1}(X\beta))\right)-b(\theta(g^{-1}(X\beta)))}{a(\phi)}+c(y,\phi).$
(2.6)

Typically maximum likelihood along with the Newton-Raphson algorithm or Fisher
scoring is used to estimate $\beta$. Usually the dispersion parameter is also
unknown and must be estimated using $\hat{\phi}$.

Suppose the inferential goal is predicting the next outcome $Y^{n+1}$. The
usual point predictor under a MSP $M$ is

$\hat{Y}^{n+1}_{M}=\hat{\mu}_{M}=g^{-1}(X^{\prime n+1}_{M}\hat{\beta}_{M}).$
(2.7)

Henceforth, our focus in on constructing valid PI’s for this point predictor.

##  3 Candidate PI’s

In this section we define four PI’s. The first is derived in our Theorem 3.1.
The second is an improvement on this interval by incorporating an extra
optimization to ensure more rapid convergence to the nominal coverage. Both of
these are in Subsec. 3.1. In Subsec. 3.2 we give our third and fourth
intervals are based on bootstrapping approach. We will argue that our
optimized interval provides the best performance.

###  3.1 Main Result and Two PI’s

One choice for a PI uses asymptotic normality of the point predictor (2.7).
Define the statistic

$Z_{pred}=Z_{pred}(M)=\frac{\hat{Y}^{n+1}_{M}-Y^{n+1}}{\sqrt{Var({\hat{Y}^{n+1}_{M}-Y^{n+1})}}}.$
(3.1)

We have the following result giving our first PI.

######  Theorem 3.1.

Suppose $Y^{n},Y^{n+1}$ come from an exponential family distribution and let
$M$ be a consistent MSP. An asymptotically normal prediction interval for a
new outcome derived from a GLM, is $PI(M)$ given by

$g^{-1}(X^{\prime n+1}_{M}\hat{\beta}_{M})\pm
z_{1-\alpha/2}\sqrt{\left[\frac{d}{d\eta}g^{-1}(\hat{\eta}^{n+1}_{M})\right]^{2}\
X^{\prime
n+1}_{M}Var(\hat{\beta}_{M})X^{n+1}_{M}+a(\hat{\phi})_{M}V(\hat{\mu})_{M}}.$
(3.2)

Proof:

Asymptotically, for any fixed $m$,

$\hat{\beta}_{m}\sim N\left(\beta_{m},(X^{\prime}_{m}WX_{m})^{-1}\right)$

where $W=(DVD)^{-1}$. In this notation, $V=\hbox{diag}[Var(y_{i})]$ is the
$n\times n$ variance matrix of observations,
$D=\hbox{diag}[\frac{\partial\eta_{i}}{\partial\mu_{i}}]$ is the $n\times n$
matrix of derivatives and $\mu$ is the $n\times 1$ mean vector. This implies

$\sqrt{n}(\hat{\beta}_{m}-\beta_{m})\overset{D}{\rightarrow}N(0,(X^{\prime}_{m_{T}}WX_{m_{T}})^{-1})$
(3.3)

where $0\in\mathbb{R}^{p}$ and
$(X^{\prime}_{m}WX_{m})^{-1}\in\mathbb{R}^{|M|\times|M|}$.

While this is useful, it is only a step toward the convergence
$\hat{\beta}_{M}\rightarrow\beta_{m_{T}}$. Since $M$ is consistent,
$m_{T}\in\mathcal{M}$ and hence $\hat{\delta_{j}}\rightarrow\delta_{j,T}$ with
probability 1 for all $j$ which implies $M\rightarrow m_{T}$. Hence, with this
assumption we get an analog to (3.3)

$\sqrt{n}(\hat{\beta}_{M}-\beta_{m_{T}})\overset{D}{\rightarrow}N(0,V^{*}_{m_{T}})$
(3.4)

where $V^{*}_{m_{T}}=(X^{\prime}_{m_{T}}WX_{m_{T}})^{-1}$. Now we define the
set

$S_{n}=\\{\omega|\forall j,\hat{\delta}_{j}(\omega)=\delta_{j,T}\\}$

and let $\textbf{1}_{S_{n}}$ be the indicator that $\omega\in S_{n}$. Further,
let $\textbf{1}_{S^{c}_{n}}$ be the indicator that $\omega$ is in the
complement of $S_{n}$ and write

$\sqrt{n}(X^{\prime n+1}_{M}\hat{\beta}_{M}-X^{\prime
n+1}_{m_{T}}\beta_{m_{T}})=\sqrt{n}(X^{\prime
n+1}_{M}\hat{\beta}_{M}-X^{\prime
n+1}_{m_{T}}\beta_{m_{T}})\textbf{1}_{S_{n}}+\sqrt{n}(X^{\prime
n+1}_{M}\hat{\beta}_{M}-X^{\prime
n+1}_{m_{T}}\beta_{m_{T}})\textbf{1}_{S^{c}_{n}}.$ (3.5)

First, note that the first term on the RHS of (3.5) becomes

$\sqrt{n}(X^{\prime n+1}_{M}\hat{\beta}_{M}-X^{\prime
n+1}_{m_{T}}\beta_{m_{T}})\textbf{1}_{S_{n}}=\sqrt{n}(X^{\prime
n+1}_{m_{T}}\hat{\beta}_{m_{T}}-X^{\prime
n+1}_{m_{T}}\beta_{m_{T}})\textbf{1}_{S_{n}}$

under consistent model selection. This term clearly converges in distribution
to a normal. Namely

$\sqrt{n}(X^{\prime n+1}_{m_{T}}\hat{\beta}_{m_{T}}-X^{\prime
n+1}_{m_{T}}\beta_{m_{T}})\textbf{1}_{S_{n}}\overset{D}{\rightarrow}N\left(0,X^{\prime
n+1}_{m_{T}}V^{*}_{m_{T}}X^{n+1}_{m_{T}}\right)$

Now observe the second term on the RHS of (3.5) can be written as

$\displaystyle\sqrt{n}(X^{\prime n+1}_{M}\hat{\beta}_{M}-X^{\prime
n+1}_{m_{T}}\beta_{m_{T}})\textbf{1}_{S^{c}_{n}}$
$\displaystyle=\sqrt{n}(X^{\prime n+1}_{M}\hat{\beta}_{M}-X^{\prime
n+1}_{M}\beta_{m_{T}}+X^{\prime n+1}_{M}\beta_{m_{T}}-X^{\prime
n+1}_{m_{T}}\beta_{m_{T}})\textbf{1}_{S^{c}_{n}}$
$\displaystyle=\sqrt{n}(X^{\prime n+1}_{M}\hat{\beta}_{M}-X^{\prime
n+1}_{M}\beta_{m_{T}})\textbf{1}_{S^{c}_{n}}+\sqrt{n}(X^{\prime
n+1}_{M}\beta_{m_{T}}-X^{\prime
n+1}_{m_{T}}\beta_{m_{T}})\textbf{1}_{S^{c}_{n}}$ $\displaystyle=X^{\prime
n+1}_{M}\sqrt{n}(\hat{\beta}_{M}-\beta_{m_{T}})\textbf{1}_{S^{c}_{n}}+\sqrt{n}(X^{\prime
n+1}_{M}-X^{\prime n+1}_{m_{T}})\beta_{m_{T}}\textbf{1}_{S^{c}_{n}}.$ (3.6)

Using (3.4) and the fact that $X^{\prime n+1}_{M}$ is bounded, we know the
first in (3.1) converges in distribution to a normal. Also, we see that in the
second term in (3.1), $\beta_{m_{T}}$ is a bounded constant vector,
$(X^{\prime n+1}_{M}-X^{\prime n+1}_{m_{T}})$ is bounded and
$P(S^{c}_{n})\overset{p}{\rightarrow}0$ by assumption. Thus we see
asymptotically that (3.5) is

$\displaystyle\sqrt{n}(X^{\prime n+1}_{M}\hat{\beta}_{M}-X^{\prime
n+1}_{m_{T}}\beta_{m_{T}})$ $\displaystyle\overset{\infty}{=}N(0,X^{\prime
n+1}_{m_{T}}V^{*}_{m_{T}}X^{n+1}_{m_{T}})\textbf{1}_{S_{n}}$
$\displaystyle+N(0,X^{\prime
n+1}_{m_{T}}V^{*}_{m_{T}}X^{n+1}_{m_{T}})\textbf{1}_{S^{c}_{n}}+\sqrt{n}\textbf{1}_{S^{c}_{n}}$

Now we need to show that $\sqrt{n}\textbf{1}_{S^{c}_{n}}=o_{p}(1)$. First note
that the union of events bound gives

$P(S^{c}_{n})\leq\sum^{p}_{j=1}P(|\hat{\delta}_{j}-\delta_{j,T}|>\eta),$

for some $\eta>0$, so using symmetry in the MSP it is enough to show

$\lim_{n\rightarrow\infty}P(|\hat{\delta}_{j}-\delta_{j,T}|>1/\sqrt{n})=0$

for any $j$. It is easy to see that

$\lim_{n\rightarrow\infty}P(|\hat{\delta}_{j}-\delta_{j,T}|>1/\sqrt{n})=\lim_{n\rightarrow\infty}P(|\hat{\delta}_{j}-\delta_{j,T}|=1)$
(3.7)

because $\delta_{j,T}$ and $\hat{\delta}_{j}$ are either 1 or 0. Now because
we have chosen a consistent MSP we have

$\lim_{n\rightarrow\infty}P(|\hat{\delta}_{j}-\delta_{j,T}|=1)=0.$

Hence, the left hand side of (3.7) is also equal to zero, implying that
$\sqrt{n}\textbf{1}_{S^{c}_{n}}=o_{p}(1)$. Now, Slutsky’s theorem gives us,

$\sqrt{n}(X^{\prime n+1}_{M}\hat{\beta}_{M}-X^{\prime
n+1}_{m_{T}}\beta_{m_{T}})\overset{D}{\rightarrow}N\left(0,X^{\prime
n+1}_{m_{T}}V^{*}_{m_{T}}X^{n+1}_{m_{T}}\right).$

Now to get a predictive distribution, we observe that the delta method gives
us

$\sqrt{n}\left(g^{-1}(X^{\prime n+1}_{M}\hat{\beta}_{M})-g^{-1}(X^{\prime
n+1}_{m_{T}}\beta_{m_{T}})\right)\overset{D}{\rightarrow}N\left(0,\left[\frac{d}{d\eta}g^{-1}(\eta^{n+1}_{m_{T}})\right]^{2}X^{\prime
n+1}_{m_{T}}V^{*}_{m_{T}}X^{n+1}_{m_{T}}\right)$ (3.8)

where $\eta^{n+1}_{m_{T}}=X^{\prime n+1}_{m_{T}}\beta_{m_{T}}$ Thus, we see
that although the GLM estimates are biased we still get convergence in
distribution when model selection occurs in the $\mathcal{M}$-closed case. The
variance of $\hat{Y}^{n+1}_{M}-Y^{n+1}$ is

$\displaystyle Var(\hat{Y}^{n+1}_{M}-Y^{n+1})$
$\displaystyle=Var(\hat{Y}^{n+1}_{M})+Var(Y^{n+1})$
$\displaystyle=Var(g^{-1}(X^{\prime n+1}_{M}\hat{\beta}_{M}))+a(\phi)V(\mu)$
$\displaystyle=\frac{1}{n}\left[\frac{d}{d\eta}g^{-1}(\eta^{n+1}_{m_{T}})\right]^{2}X^{\prime
n+1}_{m_{T}}V^{*}_{m_{T}}X^{n+1}_{m_{T}}+a(\phi)V(\mu)$ (3.9)

due to (3.8). Again, because $Y^{n+1}$ is a random variable, and not a
parameter, we must consider the variance of it as well, which we get assuming
it will come from the exponential family distribution as $Y_{1},\ldots,Y_{n}$.
This quantity, however, is impossible to compute because we do not know
$m_{T}$. Hence, we must replace $m_{T}$ with $M$, making the variance a random
quantity that depends on model selection. Now we use (3.1) as a pivotal
quantity to get

$\displaystyle 1-\alpha$ $\displaystyle\leq
P\left(\left|Z_{pred}\right|<z_{1-\alpha/2}\right)$
$\displaystyle=P\left(\left|\hat{Y}^{n+1}-Y^{n+1}\right|<z_{1-\alpha/2}\sqrt{Var(\hat{Y}^{n+1}-Y^{n+1})}\right)$
$\displaystyle=P\left(\hat{Y}^{n+1}-z_{1-\alpha/2}\sqrt{Var(\hat{Y}^{n+1}-Y^{n+1})}<Y_{n+1}<\hat{Y}^{n+1}+z_{1-\alpha/2}\sqrt{Var(\hat{Y}^{n+1}-Y^{n+1})}\right).$
(3.10)

Hence using (3.1)

$\left[g^{-1}(X^{\prime n+1}_{M}\hat{\beta}_{M})\pm
z_{1-\alpha/2}\sqrt{\frac{1}{n}\left[\frac{d}{d\eta}\widehat{g^{-1}(\eta^{n+1}_{M})}\right]^{2}\
X^{\prime
n+1}_{M}V^{*}_{M}X^{n+1}_{M}+a(\hat{\phi}_{M})V(\hat{\mu}_{M})}\right]$

is a $100(1-\alpha)\%$ prediction interval for $Y^{n+1}$

We now offer an improvement on the PI from Theorem 3.1. Note that thisn result
uses the standard normal quantile to define the predictive interval, but this
is not the only choice. Instead, we can adjust the width of the interval to
correct to poor coverage. To do this, we write the interval as

$\displaystyle PI(C_{\alpha,M})=g^{-1}(X^{\prime n+1}_{M}\hat{\beta}_{M})\pm$
$\displaystyle
C_{\alpha,M}\sqrt{\left[\frac{d}{d\eta}g^{-1}(\hat{\eta}^{n+1}_{M})\right]^{2}\
X^{\prime n+1}_{M}Var(\hat{\beta}_{M})X^{n+1}_{M}+a(\hat{\phi})V(\hat{\mu})}$
(3.11)

where $C_{\alpha,M}$ is chosen to satisfy

$C_{\alpha,M}=\arg\min_{C}P\left(Y^{n+1}\in PI(M,C)\right)$ (3.12)

for all $C$ such that $P\left(Y^{n+1}\in PI(M,C)\right)\geq 1-\alpha$.
Importantly, this probability also sees the random variable $M$ and hence
inherits the uncertainty associated with $M$ as well and $Y^{n+1}$. This is in
the same spirit of the PoSI constant in Berk et al (2013). That is, we enlarge
$C_{\alpha,M}$ to account for the uncertainty in $M$. We can approximate
$C_{\alpha,M}$ using Monte Carlo cross validation.

We begin by choosing an interval on $\mathbb{R}^{+}$ denoted
$\mathcal{C}=[C_{low},C_{high}]$ that we will perform the line search on to
estimate $C_{\alpha,M}$. Next, we randomly splitting $\mathcal{D}_{n}$ into
$L\in\mathbb{N}$ test and train sets, $\mathcal{D}_{train,\ell}$ and
$\mathcal{D}_{test,\ell}$ for $\ell=1,\ldots,L$. Then for each $\ell$ we
estimate $\beta$ by

$\hat{\beta}_{\ell}=(X^{\prime}_{train,M,\ell}X_{train,M,\ell})^{-1}X^{\prime}_{train,M,\ell}y_{train,\ell}$

using $\mathcal{D}_{train,m}$, and form the predictor
$\hat{Y}^{test}_{M,\ell}=g^{-1}(X_{test,M,\ell}\hat{\beta}_{\ell})$. Now we
form the prediction interval $PI_{M,\ell}(C)$ namely
$\hat{Y}^{n+1}_{M,\ell}\pm$

$\displaystyle
C\sqrt{\left[\frac{d}{d\eta}g^{-1}(\hat{\eta}^{n+1}_{test,M,\ell})\right]^{2}\
X^{\prime
n+1}_{test,M,\ell}Var(\hat{\beta})_{test,M,\ell}X^{n+1}_{test,M,\ell}+a(\hat{\phi})_{M,\ell}V(\hat{\mu})_{M,\ell}},$
(3.13)

and for each $C\in\mathcal{C}$, check if $y_{test,\ell}\in PI_{M,\ell}(C)$.
Then we choose the value $C$ that gives us $1-\alpha$ coverage for the Monte
Carlo samples. More formally, we can approximate $C_{\alpha,M}$ by

$\hat{C}^{MC}=\arg\min_{C}\frac{1}{L}\sum^{L}_{\ell=1}\left|\frac{1}{\\#(\mathcal{D}_{test,\ell})}\sum^{\\#(\mathcal{D}_{test,\ell})}_{i=1}I_{y_{test,\ell}\in
PI_{M,\ell}(C)}-(1-\alpha)\right|$ (3.14)

where $I_{y_{test,\ell}\in PI_{M,\ell}(C)}$ is the indicator that the test
values are in the constructed intervals.

The intuition behind using this interval in place of the PI in Theorem 3.1 is
that, in finite samples the difference between $z_{1-\alpha/2}$ and
$\hat{C}^{MC}$ can be interpreted as t he added variability due to model
uncertainty.

###  3.2 Two Bootstrap Based PI’s

In the frequentist setting, perhaps the most natural way to obtain a
prediction interval that takes into account the uncertainty of both model
selection and the uncertainty associated with the distribution of the new
outcome is to make use of the bootstrap. Accordingly, to form our fiest
bootstraped PI, we use the bootstrap to estimate the distribution of

$\hat{\mu}_{M}=E(Y^{n+1}|X^{n+1}_{M})=g^{-1}(X^{n+1}_{M}\hat{\beta}_{M}),$
(3.15)

and $a(\hat{\phi})_{M}$. Then for each bootstrapped mean and dispersion
function, we generate a new observation from the distribution of
$Y^{n+1}|X^{n+1},\mu,\phi$, i.e. $\mathcal{G}$.

Let $\hat{p}(\hat{\mu})$ denote be the bootstrapped density of (3.15) and
$\hat{p}(a(\hat{\phi}))$ be the bootstrapped density of $a(\hat{\phi})_{M}$.
Then $\hat{p}(Y^{n+1})$ be the resulting estimated density of $Y^{n+1}$.

The procedure is as follows,

  * •

obtain $B$ bootstrap replications of $\hat{\mu}_{M}$, denoted
$\mu^{*}_{1},\ldots,\mu^{*}_{B}$

  * •

obtain $B$ bootstrap replications of $a(\hat{\phi})_{M}$, denoted
$a(\phi)^{*}_{1},\ldots,a(\phi)^{*}_{B}$

  * •

generate
$y^{*}_{1}(\mu^{*}_{1},a(\phi)^{*}_{1}),\ldots,y^{*}_{B}(\mu^{*}_{B},a(\phi)^{*}_{B})$,
from $\mathcal{G}$.

The sample $y^{*}_{1},\ldots,y^{*}_{B}$ can be used to estimate an approximate
marginal predictive distribution for $Y^{n+1}$. To obtain the PI, we use the
appropriate percentile interval from this distribution. That is, to obtain a
$100(1-\alpha)\%$ PI we use the interval

$[q^{*}_{1-\alpha/2},q^{*}_{\alpha/2}]$ (3.16)

where $q^{*}_{\alpha}$ is the $\alpha$ quantile from $\hat{p}(Y^{n+1})$. The
use of $\hat{p}(\hat{\mu})$ and $\hat{p}(a(\hat{\phi}))$ to obtain the
estimated predictive distribution $\hat{p}(Y^{n+1})$ allows $\hat{p}(Y^{n+1})$
to inherit the variability from $\hat{p}(\hat{\mu})$, $\hat{p}(a(\hat{\phi}))$
and the variability that is already associated with the known parametric
distribution $\mathcal{G}$. Hence, the interval (3.16) is typically widened
due to the uncertainty of the model selection procedure as well as the
uncertainty of the distribution of $Y^{n+1}$.

Now, in the GLM setting, coverage for the PI in (3.16) should be closer to the
$1-\alpha$ nominal coverage than the PI resulting from ignoring the model
uncertainty. Note that as $n\rightarrow\infty$ the variability due to model
uncertainty will go to 0 and this interval will converge to the standard PI.
Bootstrap PI’s for the Gaussian case are studied in a fairly narrow (small $d$
and moderate $n$) setting in Hong et al (2018b). These authors suggest that in
this setting the bootstrap distribution fails to assess the uncertainty of
model selection accurately. We explore different simulation settings to
evaluate the performance of bootstrap intervals in Subsec. 4.

Our second bootstrapped PI is formed as follows. Recall he interval in Theorem
3.1 is a random because it depends on $M$. It is directly usable for
predictions, but we must use $\hat{M}$ in place of $M$ to get a confidence
statement. Nevertheless, we provide an approximate interval by “smoothing”
over $M$, which accounts for the uncertainty of $M$ in both the center and
width of the interval. This is similar to the approach used in Efron (2014)
for estimation. The method we propose is to use $\hat{p}(\hat{\mu})$, the
bootstrap the distribution of $\hat{\mu}_{M}=g^{-1}(\hat{\eta}_{M})$ as
described earlier in this Subsection, to obtain an approximation for the
predictor and its variance that accounts for model selection uncertainty.

Specifically, we use

$\tilde{\mu}=\frac{1}{B}\sum^{B}_{b=1}\mu^{*}_{b}$

for the point predictor. We approximate the variance of $\hat{\mu}_{M}$ with

$\displaystyle Var(\mu^{*})$ $\displaystyle=\widehat{Var}(\hat{\mu}_{M})$
$\displaystyle=\frac{1}{B-1}\sum^{B}_{b=1}\left(\mu^{*}_{b}-\tilde{\mu}\right)^{2}$
$\displaystyle\approx\left[\frac{d}{d\eta}g^{-1}(\hat{\eta}^{n+1}_{M})\right]^{2}\
X^{\prime n+1}_{M}Var(\hat{\beta}_{M})X^{n+1}_{M},$

and the estimated variance of the predictive distribution is given by

$\displaystyle Var(Y^{*})$ $\displaystyle=\widehat{Var}(Y^{n+1})$
$\displaystyle=\frac{1}{B-1}\sum^{B}_{b=1}\left(y^{*}(\mu^{*}_{b})-\bar{y}^{*}\right)^{2}$
$\displaystyle\approx a(\hat{\phi}_{M})V(\hat{\mu}_{M})$

where $\bar{y}^{*}=\frac{1}{B-1}\sum^{B}_{b=1}y^{*}(\mu^{*}_{b})$. We treat
$Y^{*}$ as a random variable approximating $Y^{n+1}$. Note also that we are
required to estimate $a(\hat{\phi}_{M})V(\hat{\mu}_{M})$ in (3.2), but this
again is a random variable so we using bootstrapping to account for the
uncertainty in $M$ is necessary for this term also. Now as an ad-hoc fix, we
rewrite (3.2) to give our second bootstrapped PI

$PI(M)=\tilde{\mu}\pm z_{1-\alpha/2}\sqrt{Var(\mu^{*})+Var(Y^{*}))}.$ (3.17)

##  4 Simulation Results for GLM’s

We give two contexts in which the PI’s we have defined in (3.2), (3.11),
(3.16), (3.17) can be readily found. Respectively, these intervals are labeled
the asymptotic normal PI (AN), the optimized AN $\hat{C}^{MC}$, the
bootstrapped (boot) PI, and the ‘smoothed’ asymptotic normal (S-AN) PI. In
addition to the intervals we have derived, we give the BMA PI’s as well as the
‘Naive’ PI’s obtained by applying the inverse link to a confidence interval
for the mean on the linear predictor scale; this is often done by
practitioners as a pragmatic solution.

In Sec. 4.1 we present these intervals for the standard Gaussian case and in
Sec. 4.2, we present the prediction intervals for binomial regression, i.e., a
more general case of logistic regression. For both cases we use 500 new
observations from their respective distribution and calculate the estimated
predictive coverage using

$\widehat{coverage}=\frac{1}{500}\sum^{500}_{i=1}I_{y^{new}_{i}\in
PI_{i}(X^{new}_{i},X^{n},y^{n})},$ (4.1)

where each $PI_{i}(X^{new}_{i},X^{n},y^{n})$ depends on the data and the new
observed explanatory variables. For the PIs that require bootstrapping we
resample the data 500 times to obtain the bootstrapped distributions.

###  4.1 Gaussian Linear Models

In the standard case, we assume $Y\sim N(\mu,\sigma^{2})$, and the log
likelihood is

$L(\mu_{i},\sigma^{2}|y_{i})=\frac{y_{i}\mu_{i}-(\mu_{i}^{2}/2)}{\sigma^{2}}-\left(\frac{y^{2}_{i}}{2\sigma^{2}}+\log(\sigma\sqrt{2\pi}\right),$

the canonical parameter is $\theta_{i}=\mu_{i}$, $b(\theta)=\mu+i^{2}/2$,
$a(\phi)=\sigma^{2}$, and $V(\mu_{i})=1$. The linear predictor uses the
identity link function and the point predictor is

$\hat{Y}^{n+1}_{M}=X^{\prime n+1}_{M}\hat{\beta}_{M}.$

The asymptotic normal PI from (3.2) for $Y^{n+1}$ is

$PI(M)=\left[\hat{Y}^{n+1}_{M}\pm
z_{1-\alpha/2}\hat{\sigma}_{M}\sqrt{X^{\prime
n+1}_{M}(X^{\prime}_{M}X_{M})^{-1}X^{n+1}_{M}+1}\right].$

Our simulation results for Gaussian data includes coverage and width estimates
for the normal PI in (3.2), the PI (3.13) using $\hat{C}^{MC}$, the bootstrap
PI in (3.16), and the ‘smoothed’ normal interval (3.17). We do not include the
Naive interval because in the Gaussian case it is equivalent to AN. For the
interval using $\hat{C}^{MC}$, we do a grid search for the value of
$C_{\alpha,M}$ on the interval from 1.95 to 5 in increments of 0.05.

The simulation setup is as follows. First, we consider two model selection
procedures, BIC and AIC. Both methods are implemented in R using the step()
function by setting the respective penalties for BIC and AIC. We also use BMA
implemented with the BAS package in R. We consider various choices for $n$
(30,50,100,200) and choose $p=25$. We randomly generate values for $\sigma$
and $\beta$ once, and fix those values throughout the simulations.
Accordingly, let

$\beta=(\beta_{1},\ldots,\beta_{25})^{\prime}=(6.43,4.39,4.26,4.11,0,\ldots,0)^{\prime}$

and $\sigma=0.93$. We simulate $n$ observations for the design matrix $X$
according to

$X\sim MVN_{p}(0,I_{p}),$

and then draw and $n\times 1$ vector of observations from $Y\sim
N(X\beta,\sigma^{2}I_{n})$. We then calculate estimated coverage using (4.1).
Ideally, we want coverage close to $0.95$. When choosing between competing
PI’s with good coverage, we prefer the one with the narrowest width. The
results are seen in Table 1.

Table 1: Simulation results for Gaussian data with $p=25$ and $p_{0}=4$.  $n$ | MSP | Interval | Coverage | Avg.Width (SE)
---|---|---|---|---
50 | AIC | AN | 0.88 | 3.63 (0.17)
|  | S-AN | 1 | 14.2 (3.81)
|  | boot | 0.99 | 8.8 (1.87)
|  | $\hat{C}^{MC}$ | 0.98 | 5.5 (0.35)
| BIC | AN | 0.88 | 3.63 (0.18)
|  | S-AN | 1 | 13.79 (3.80)
|  | boot | 0.99 | 8.4 (1.78)
|  | $\hat{C}^{MC}$ | 0.98 | 5.37 (0.26)
| BMA |  | 0.89 | 3.98 (0.11)
100 | AIC | AN | 0.91 | 3.63 (0.08)
|  | S-AN | 0.95 | 4.43 (0.41)
|  | boot | 0.92 | 3.81 (0.26)
|  | $\hat{C}^{MC}$ | 0.93 | 4.08 (0.09)
| BIC | AN | 0.92 | 3.70 (0.06)
|  | S-AN | 0.94 | 4.25 (0.33)
|  | boot | 0.92 | 3.74 (0.22)
|  | $\hat{C}^{MC}$ | 0.94 | 3.97 (0.05)
| BMA |  | 0.93 | 3.73 (0.05)
200 | AIC | AN | 0.92 | 3.70 (0.05)
|  | S-AN | 0.93 | 3.95 (0.16)
|  | boot | 0.91 | 3.64 (0.14)
|  | $\hat{C}^{MC}$ | 0.93 | 3.96 (0.06)
| BIC | AN | 0.94 | 3.75 (0.03)
|  | S-AN | 0.94 | 3.91 (0.14)
|  | boot | 0.92 | 3.65 (0.12)
|  | $\hat{C}^{MC}$ | 0.94 | 3.92 (0.03)
| BMA |  | 0.92 | 3.74 (0.03)
Table 2: Gaussian cross validation results for the optimal choice for $\hat{C}^{MC}$. n | MSP | $\hat{C}^{MC}$
---|---|---
50 | AIC | 2.95
| BIC | 2.90
100 | AIC | 2.20
| BIC | 2.10
200 | AIC | 2.10
| BIC | 2.05

Note that the differences between using AIC and BIC are negligible, so we
describe the performance of each PI only once (rather than once for each MSP).
It is seen in Table 1 that AN has low coverage for $n=50$, but gets close to
the nominal coverage for the larger sample sizes. For $n=50$, both S-AN and
boot give at least the nominal coverage and arguably reasonable width of PI’s
to be useful. Here, $\hat{C}^{MC}$ gives close to the stated 95% coverage and
is noticeably narrower than both S-AN and boot, so it is the preferred PI.

When $n=100$ and 200, we observe all of the 5 PIs are roughly equal in terms
of coverage and width. Since AN is the easiest to implement as it does not
require any bootstrapping or cross validation, we recommend using it with
relatively large $n$. For intermediate $n$ we recommend using $\hat{C}^{MC}$
as it gives appropriate coverage and is narrower than the other PIs.

We give the optimal choices for $\hat{C}^{MC}$ for each sample size in Table
2. We observe that as sample size increases, $\hat{C}^{MC}$ decreases as
expected. This reflects the fact that as we gather more data, the uncertainty
in model selection also decreases.

###  4.2 Binomial Regression

Suppose we have $n$ independent but not identically distributed random
variables following $Y_{i}\sim Bin(r_{i},p_{i})$ so $E(Y_{i})=r_{i}p_{i}$. We
write $W=\frac{Y_{i}}{r_{i}}$ as our response to model the proportion of
success, and then we convert back to number of successes to form our
predictive interval. Now we have $E(W)=p_{i}$ and the log likelihood for a
given $i$ is given by

$L(p_{i}|w_{i})=\frac{w_{i}\log\left(\frac{p_{i}}{1-p_{i}}\right)+\log(1-p_{i})}{\frac{1}{r_{i}}}+\log\binom{r_{i}}{nw_{i}},$

which reveals the canonical parameter

$\theta_{i}=logit(p_{i})=\log\left(\frac{p_{i}}{1-p_{i}}\right).$

We also see that $a(\phi)=\frac{1}{r_{i}}$,
$b(\theta_{i})=-\log(1+e^{\theta_{i}})=-\log(1-p_{i})$, and thus

$V(p_{i})=\frac{\partial^{2}b(\theta_{i})}{\partial
p_{i}^{2}}=\frac{p_{i}(1-p_{i})}{r_{i}}.$

Thus the linear predictor is defined by the logit link as

$E\left(\frac{Y_{i}}{r_{i}}\right)=g(p_{i})=\log\left(\frac{p_{i}}{1-p_{i}}\right)=X^{\prime}_{i}\beta$

and the inverse link function, which gives the probability of success, is
given by

$p_{i}=g^{-1}(X^{\prime}_{i}\beta)=\frac{1}{1+e^{-X^{\prime}_{i}\beta}}.$

Of course, we do not know $p_{i}$, so we estimate $p_{i}$ by

$\hat{p}_{i}=g^{-1}(X^{\prime}_{i}\hat{\beta})=\frac{1}{1+e^{-X^{\prime}_{i}\hat{\beta}}}.$

Given $n$ observations $Y_{1},\ldots,Y_{n}$, our goal is to predict the total
number of successes $Y_{n+1}$ in $r_{n+1}$ trials while accounting for model
selection. We denote the predicted probability of success $\hat{p}^{n+1}_{M}$
and its value is given by

$\hat{p}^{n+1}_{M}=\frac{1}{1+e^{-X^{\prime n+1}_{M}\hat{\beta}_{M}}}.$

Recalling that

$E(Y_{n+1})=r_{n+1}\cdot g^{-1}(X^{\prime n+1}\beta)=r_{n+1}\cdot p_{n+1},$

the form of the post-model selection AN PI for a binomial random variable is

$\displaystyle PI(M)=r_{n+1}\cdot\hat{p}^{n+1}_{M}\pm$ $\displaystyle
z_{1-\alpha/2}\cdot
r_{n+1}\sqrt{\frac{e^{-2\hat{\eta}^{n+1}_{M}}}{\left(1+e^{-\hat{\eta}^{n+1}_{M}}\right)^{4}}X^{\prime
n+1}_{M}Var(\hat{\beta}_{M})X^{n+1}_{M}+\frac{1}{r^{n+1}}\hat{p}_{M}\left(1-\hat{p}_{M}\right)}$
(4.2)

where the factor $r_{n+1}$ in the width of the intervals comes from the the
distribution in (3.8) being multiplied by this factor. The interval in (4.2)
gives a prediction interval for total number of successes in $r_{n+1}$ trials.

In the setting described above, our simulations are as follows. Let $X\sim
MVN_{p}(0,I_{p})$ and

$\beta=(\beta_{1},\ldots,\beta_{25})^{\prime}=(0.252,0.171,-0.268,0.09,0,\ldots
0)^{\prime}.$

Now we calculate the estimated coverage using (4.1). Again, we want coverage
close to $0.95$ and narrow width. The simulated results are given in Table 3.

Table 3: Simulation results for binomial data with $r=30$.  $n$ | MSP | Interval | Coverage | Avg Width (SE)
---|---|---|---|---
50 | AIC | Naive | .51 | 4.9 (1.35)
|  | AN | 0.83 | 9.51 (0.93)
|  | S-AN | 1 | 28.31 (2.22)
|  | boot | 1 | 20.82 (3.35)
|  | $\hat{C}^{MC}$ | 0.97 | 15.55 (0.78)
| BIC | Naive | .48 | 4.03 (1.08)
|  | AN | 0.90 | 9.56 (0.76)
|  | S-AN | 1 | 23.36 (3.93)
|  | boot | 1 | 17.81 (3.06)
|  | $\hat{C}^{MC}$ | 0.98 | 14.35 (0.74)
100 | AIC | Naive | .46 | 3.29 (0.86)
|  | AN | 0.94 | 9.42 (0.87)
|  | S-AN | 0.99 | 14.11(2.22)
|  | boot | 0.99 | 11.99(1.66)
|  | $\hat{C}^{MC}$ | 0.99 | 13.01 (0.61)
| BIC | Naive | .43 | 2.77 (0.69)
|  | AN | 0.95 | 9.38 (0.79)
|  | S-AN | 1 | 13.38 (2.00)
|  | boot | 0.99 | 11.57 (1.47)
|  | $\hat{C}^{MC}$ | 0.99 | 12.67 (0.54)
200 | AIC | Naive | .37 | 2.43 (0.90)
|  | AN | 0.94 | 8.91 (1.60)
|  | S-AN | 0.97 | 10.08(2.10)
|  | boot | 0.97 | 9.37 (1.96)
|  | $\hat{C}^{MC}$ | 0.94 | 9.15 (0.87)
| BIC | Naive | .37 | 2.18 (0.79)
|  | AN | 0.94 | 8.93 (1.56)
|  | S-AN | 0.97 | 9.75 (1.92)
|  | boot | 0.98 | 9.30 (1.82)
|  | $\hat{C}^{MC}$ | 0.95 | 9.01 (0.83)
Table 4: Binomial cross validation results for the optimal choice for $\hat{C}^{MC}$. n | MSP | $\hat{C}^{MC}$
---|---|---
50 | AIC | 5.00
| BIC | 4.40
100 | AIC | 4.10
| BIC | 3.90
200 | AIC | 3.00
| BIC | 2.85

For $n=50$, the Naive interval has very poor coverage for both AIC and BIC.
Using AIC and the AN interval results in undercoverage, but it is much better
than the Naive PI. This is also true using BIC as the MSP. Both S-AN and boot
are conservative, give coverage larger than the stated coverage. The width of
both S-AN and boot make the intervals fairly uninformative despite having
better coverage than Naive and AN. Finally, we observe $\hat{C}^{MC}$ performs
noticeably better than the other PI’s. This suggests that the cross validation
step to widen the asymptotic normal PI is useful.

Looking at the $n=100$ and case, we see the Naive interval is worse than the
smaller sample case. AN gives very good coverage and the smallest width among
all of the PIs for both AIC and BIC. The other 3 PIs, S-AN, boot, and
$\hat{C}^{MC}$ give close stated coverage but they are slightly wider than the
AN interval. When $n=200$ we see Naive is by far the worst among the 5 PIs,
but the other 4 are roughly the same with AN and $\hat{C}^{MC}$ having perhaps
slightly better coverage and narrower PIs than S-AN and boot.

These results confirm two main points. First, the AN PI achieves the stated
95% coverage as given in Theorem 3.1 when the sample size is large enough.
Second $\hat{C}^{MC}$ always gives appropriate coverage, and appears to reduce
to AN as $n$ increases. This leads us to recommend using $\hat{C}^{MC}$,
especially for intermediate sample sizes, and use AN for large $n$.

Again, we list the optimal cross validation constants in Table 4. As in the
Gaussian case, we see that $\hat{C}^{MC}$ decreases as the sample size
increases. However, in the Binomial case, $\hat{C}^{MC}$ is noticeably larger
than the Gaussian case. This may be due to the fact that we are using a normal
PI for data that is not normal.

##  5 Discussion

Our main contribution is the PI in Theorem 3.1, and the small sample
correction using $\hat{C}^{MC}$ given in (3.13). Much of the literature on GLM
prediction has focused on confidence intervals around predictors, which we
refer to as the ‘Naive’ PI, rather than true prediction intervals. That is, it
is common for analysts to apply the inverse link function to the endpoints of
a confidence interval in the linear predictor scale. This approach does not
account for uncertainty appropriately because it uses the variability on the
linear predictor scale rather than the data scale. Here we have presented
prediction intervals that are derived on the model-scale rather than the
linear-predictor scale.

We have presented several prediction intervals that consider model
uncertainty. The PI derived in Theorem 3.1 accounts for model uncertainty via
the consistency of the MSP. This PI severely underperforms in terms of
predictive coverage in small sample size, e.g. $n\approx p$, cases but as
$n\rightarrow\infty$ the predictive coverage is roughy the nominal $1-\alpha$
coverage. The boot and S-AN PIs tend to be too wide, suggesting far too much
model uncertainty, for these PIs to be useful. These two PIs overcorrect the
width of the intervals for the amount of the uncertainty in the model
selection considered here. Again as $n$ increases, both boot and S-AN become
usable (due to the MSP choosing the correct model). At this point, the
bootstrapping is not necessary, however. Since AN performs well with large
samples, there is no need to bootstrap, we can directly use AN.

Taken together, our results provide valid post model selection PIs for GLM’s
for moderate and large samples.

##  Appendix A Extension to GLMM’s

The approach described in Sec. 2.2 to obtain valid prediction intervals after
model selection extends naturally to the class of generalized linear mixed
models. Here we assume the random variable $Y|X,\beta,Z,U\sim{\mathcal{G}}$
where ${\mathcal{G}}$ is a distribution in the exponential family. We write
the linear predictor as

$\eta=g(E(Y|U))=g(\mu|U)=X\beta+ZU$ (5.1)

where $\beta$ is the vector of fixed effects, $U$ is a random effect such that
$U\sim N(0,\Sigma_{U})$. $X$ and $Z$ are their respective design matrices. The
mean function is

$\mu=E(Y|U)=g^{-1}(\eta)=g^{-1}(X\beta+ZU)$

and the variance is

$Var(Y|U)=V^{1/2}_{\mu}AV^{1/2}_{\mu}$

where $V^{1/2}_{\mu}=\hbox{diag}\left[\sqrt{V(\mu)}\right]$ and
$A=\hbox{diag}\left[1/a(\phi)\right]$.

As with GLM’s, model selection is often performed when forming predictors. In
the GLMM setting model selection can be done on both $X$ and $Z$, however here
we focus on model selection on the design matrix $X$. Analogous to the GLM
case, we state an asymptotic normal predictive interval For GLMM’s that is
derived in the same way as the GLMM. The only difference is the random effects
part of the linear predictor. However, recall the random effects have
expectation 0, so the location of the asymptotic distribution does not change.
The variance, on the other hand, does increase. This is seen in (5.2 ) as the
width has an extra term for the variance of the random effects: We get that
$PI(M,C_{\alpha})=g^{-1}(\hat{\eta}^{n+1}_{M})\pm$

$C_{\alpha}\sqrt{\left[\frac{d}{d\eta}g^{-1}(\hat{\eta}^{n+1}_{M})\right]^{2}\left(X^{\prime
n+1}_{M}Var(\hat{\beta}_{M})X^{n+1}_{M}+Z^{\prime
n+1}Var(\hat{u})Z^{n+1}\right)+a(\hat{\phi})V(\hat{\mu})}.$ (5.2)

##  Appendix B GLMM bootstrap intervals

As with the GLM AN interval, we can approximate the variance of
$g^{-1}(\hat{\eta}^{n+1}_{M})$ using bootstrapping and replace (5.2) with

$PI(M,C_{\alpha})=g^{-1}(\hat{\eta}^{n+1}_{M})\pm
z_{1-\alpha/2}\sqrt{\hat{Var}(g^{-1}(\hat{\eta}^{n+1}_{M}))^{boot}+a(\hat{\phi})V(\hat{\mu})}$
(5.3)

where $\hat{Var}(g^{-1}(\hat{\eta}^{n+1}_{M}))^{boot}$ is simply the variance
of the bootstrapped distribution of $g^{-1}(\hat{\eta}^{n+1}_{M})$.

We can also use the bootstrap approach, in same way as with the GLM, to obtain
a bootstrap distribution for a new outcome. In the GLMM setting, we bootstrap
the expected value of the distribution for a new outcome,

$\hat{\mu}_{M}=g^{-1}(X^{\prime n+1}_{M}\hat{\beta}+Z^{\prime n+1}\hat{u}).$

We proceed as follows:

  * •

obtain $B$ bootstrap replications of $\hat{\mu}_{M}$, denoted
$\mu^{*}_{1},\ldots,\mu^{*}_{B}$,

  * •

obtain $B$ bootstrap replications of $a(\hat{\phi})_{M}$, denoted
$a(\phi)^{*}_{1},\ldots,a(\phi)^{*}_{B}$,

  * •

generate
$y^{*}_{1}(\mu^{*}_{1},a(\phi)^{*}_{1}),\ldots,y^{*}_{B}(\mu^{*}_{B},a(\phi)^{*}_{B})$,
from $\mathcal{G}$.

The sample $y^{*}_{1},\ldots,y^{*}_{B}$ is used to obtain the predictive
interval by extracting the appropriate percentile interval from this
distribution. Thus the PI is

$[q^{*}_{1-\alpha/2},q^{*}_{\alpha/2}],$ (5.4)

the $1-\alpha/2$ and $\alpha/2$ quantiles from $y^{*}_{1},\ldots,y^{*}_{B}$
which inherits the uncertainty of $M$, $\hat{\beta}$ and $\hat{u}$. These
intervals are implementable assuming we already have estimates $\hat{\beta}$
and $\hat{u}$. Regardless of which method is used to form predictors, we
theoretically can use both intervals (5.3) or (5.4) because predictors and
mean and variance functions, as well as the uncertainty associated with the
MSP can be obtained through the bootstrap procedure. Thus, a closed form
solution of parameter estimates is not necessary to obtain valid PI’s. Also,
we can, at least theoretically, still use both of the intervals presented to
account for the uncertainty of model selection in the random effects design
matrix.

##  Appendix C Computational issues for GLMM’s

The theoretical and bootstrap based intervals we havew proposed to capture the
uncertainty of model selection are not implementable, at least yet. This is
due to convergence issues with implementing GLMM’s. Estimation in GLMM’s
requires integrating out the random effects, and these integrals do not have
closed form solutions. Thus, numerical integration is necessary, making the
integrals computationally hard.

In practice, now, there is no a single best approach so one tries many
approaches until the algorithm converges. Once convergence is achieved,
classical approaches to assess model fit are used. In the bootstrapping
approach, we require estimation over many repeated samples of the data and
this would require convergence of the estimates in the GLMM over each
resample. The estimates require numerical integration for each resample of the
data which requires a person trying several algorithms until one works. We
attempted this, but we were unsuccessful because convergence in each resample
using a fixed numerical integration method is not feasible.

Acknowledgments

The first author acknowledges funding from the University of Nebraska Program
of Excellence in Computational Science.

## References

  *   * Bachoc et al (1988) Bachoc F, Leeb H, Pötscher B (1988) Asymptotic properties of information theoretic methods of model selection. University of California at Davis, Division of Statistics Technical Report
  * Bachoc et al (2019) Bachoc F, Leeb H, Pötscher B (2019) Valid confidence intervals for post-model-selection predictors. Annals of Statistics 47:1475–1504
  * Berk et al (2013) Berk R, Brown L, A. B, et al (2013) Valid post-selection inference. Annals of Statistics 41:802–837
  * Efron (2014) Efron B (2014) Estimation and accuracy after model selection. Journal of the American Statistical Association 109:991–1022
  * Hong et al (2018a) Hong L, Kuffner T, Martin R (2018a) On overfitting and post-selection uncertainty assessments. Biometrika 105:221–224
  * Hong et al (2018b) Hong L, Kuffner T, Martin R (2018b) On prediction of future insurance claims when the model is uncertain. Preprint
  * Kuchibhotla et al (2020) Kuchibhotla AK, Brown LD, Buja A, et al (2020) A model free perspective for linear regression: Uniform-in-model bounds for post selection inference. Annals of Statistics 49:2953–2981
  * Leeb (2009) Leeb H (2009) Conditional predictive inference post model selection. Annals of Statistics 37:2838–2876
  * Leeb et al (2015) Leeb H, Pötscher B, Ewald K (2015) On various confidence intervals post-model-selection. Statistical Science 30:216–227
  * Stine (1985) Stine R (1985) Bootstrap prediction intervals for regression. Journal of the American Statistical Association 80:1026–1031
  * Stine (2021) Stine R (2021) Prediction intervals for glms, gams, and some survival regression models. Communications in Statistics - Theory and Methods


# From Words to Routes: Applying Large Language Models to Vehicle Routing

Zhehui Huang, Guangyao Shi, Gaurav S. Sukhatme The authors are with the
Department of Computer Science, University of Southern California, Los
Angeles, CA 90089 USA (e-mail: zhehuihu, shig, gaurav@usc.edu). GSS holds
concurrent appointments as a Professor at USC and as an Amazon Scholar. This
paper describes work performed at USC and is not associated with Amazon.

###### Abstract

LLMs have shown impressive progress in robotics (e.g., manipulation and
navigation) with natural language task descriptions. The success of LLMs in
these tasks leads us to wonder: What is the ability of LLMs to solve vehicle
routing problems (VRPs) with natural language task descriptions? In this work,
we study this question in three steps. First, we construct a dataset with 21
types of single- or multi-vehicle routing problems. Second, we evaluate the
performance of LLMs across four basic prompt paradigms of text-to-code
generation, each involving different types of text input. We find that the
basic prompt paradigm, which generates code directly from natural language
task descriptions, performs the best for GPT-4, achieving $56\%$ feasibility,
$40\%$ optimality, and $53\%$ efficiency. Third, based on the observation that
LLMs may not be able to provide correct solutions at the initial attempt, we
propose a framework that enables LLMs to refine solutions through self-
reflection, including self-debugging and self-verification. With GPT-4, our
proposed framework achieves a $16\%$ increase in feasibility, a $7\%$ increase
in optimality, and a $15\%$ increase in efficiency. Moreover, we examine the
sensitivity of GPT-4 to task descriptions, specifically focusing on how its
performance changes when certain details are omitted from the task
descriptions, yet the core meaning is preserved. Our findings reveal that such
omissions lead to a notable decrease in performance: $4\%$ in feasibility,
$4\%$ in optimality, and $5\%$ in efficiency. Website:
https://sites.google.com/view/words-to-routes/

##  I Introduction

Vehicle Routing Problems (VRPs) with numerous real-world practical
applications, including transportation [1], e-commerce [2], and robotics [3],
has been discussed since last century [4]. However, despite the existence of
numerous algorithms for solving various VRPs [5, 6, 7], no existing work has
yet provided feasible, let alone optimal, solutions based entirely on Natural
Language (NL) task descriptions. Consequently, whenever a new VRP variant
arises, we rely on experts to identify and modify an appropriate existing
algorithm or need experts to develop a new algorithm to solve the problem.
Both approaches heavily depend on domain expertise and require significant
manual effort.

Recently, Large Language Models (LLMs) have shown remarkable capabilities in
solving robotic problems with NL task descriptions, such as manipulation [8,
9, 10] and navigation [11, 12, 13]. Such progress motivates us to explore
using LLMs to solve VRPs, which are widely used to model robotic decision-
making problems. Inspired by the observation that LLMs perform poorly when
directly generating solutions from NL task descriptions [14] and the
remarkable capability of LLMs in text-to-code generation [15], we adopt the
method that asks LLMs to solve VRPs by generating code given NL task
descriptions (Figure 1).

Figure 1: Overview: LLMs for vehicle routing problems. The task description is
from users, and environment information is given. Figure 2: Framework
overview. LLM $\\#1$, LLM $\\#2$, and LLM $\\#3$ are the same model but with
different roles in our settings.

In this work, we investigate the performance of pre-trained LLMs for VRPs
using zero-shot prompting and without finetuning LLMs. To this end, we first
construct a dataset with 21 variants of VRPs (10 for a single vehicle and 11
for multiple vehicles) and evaluate the performance of LLMs across all tasks
in this dataset with three metrics: feasibility, optimality, and efficiency.
In each VRP variant, we initially provide a thorough task description and aim
to reduce the ambiguity from the user’s perspective. After constructing the
dataset, we aim to identify an effective prompt paradigm to fully utilize
LLMs. Our search process has two stages. In the first stage, we propose and
evaluate four prompt paradigms to assess the task-understanding ability and
problem-solving capabilities of LLMs. We find that directly feeding NL to LLMs
achieves the best performance for GPT-4 and Gemini 1.0 Pro. GPT-4 achieves
$56\%$ feasibility, $40\%$ optimality, and $53\%$ efficiency, while Gemini 1.0
Pro only achieves $1\%$ feasibility, $0\%$ optimality, and $1\%$ efficiency.
Moreover, we observe that generated programs might not be executable or
provide the correct solution at the first attempt, which limits the
performance of LLMs.

To increase the performance of LLMs in providing feasible or even optimal
solutions, we propose a framework by asking LLMs to regenerate solutions.
Figure 2 shows an overview of our framework. This framework benefits GPT-4 to
achieve a $16\%$ increase in feasibility, a $7\%$ increase in optimality, and
a $15\%$ increase in efficiency, while it does not increase the performance of
Gemini 1.0 Pro.

In addition to proposing a framework to make LLMs more powerful, we further
study the question: What if the NL task descriptions from users are not
thorough? Such studies can help us understand the sensitivity of LLMs to the
NL task description. We build a new dataset by deleting details from NL task
descriptions while ensuring their meanings remain unchanged. The performance
of GPT-4 in this new dataset drops $4\%$ in feasibility, $4\%$ in optimality,
and $5\%$ in efficiency compared with the original dataset. However,
identifying what details are necessary for LLMs is challenging. Instead, we
propose a simple but efficient mechanism that asks LLMs to produce thorough NL
task descriptions by interacting with users.

To summarize, our contributions are as follows:

  * •

We construct a dataset that can be used as a benchmark for LLMs in VRPs.

  * •

We show that feeding NL task descriptions directly into LLMs can achieve the
best performance without any finetuning of LLMs. By only feeding NL task
descriptions into LLMs and asking LLMs to generate Python code to solve
vehicle routing problems, GPT-4 achieves $56\%$ feasibility, $40\%$
optimality, and $53\%$ efficiency in our dataset.

  * •

We propose a framework that improves the performance of GPT-4 with a $16\%$
increase in feasibility, a $7\%$ increase in optimality, and a $15\%$ increase
in our dataset. We propose a mechanism to improve the performance of LLMs when
NL task descriptions are not thorough.

##  II Related Work

###  II-A LLMs for Vehicle Routing Problems

There is very limited work that uses LLMs to solve VRPs. However, there are
some works that focus on integrating LLMs into the optimization algorithms use
the classical single traveling salesman problem (TSP) as a test case [16, 17,
18, 19] or use the capacitated VRPs as test cases [20]. Yang et al. [16] use
LLMs as optimizers and iteratively generate new solutions given task
descriptions and previous solution-score pairs. However, this paper relies on
initial solution-score pairs to construct the initial input prompt, making it
not rely solely on task descriptions. Huang et al. [20] use multimodal LLMs to
solve optimization problems and rely on XML data as additional input to solve
the problem. Liu et al. [17], Liu et al. [18], and Ye et al. [19] focus on
combining LLMs with evolutionary algorithms. Liu et al. [17] ask LLMs to
select parent solutions and perform crossover and mutation to generate
offspring solutions. Liu et al. [18], instead of asking LLMs to operate on
solutions [17], ask LLMs to operate on heuristic algorithms. Ye et al. [19]
share the same idea as [18], but it introduces short-term reflection and long-
term reflection mechanisms to use history heuristic algorithms better.
However, all these methods do not solely rely on NL task descriptions.
Besides, these methods are based on the idea that iteratively improves the
solutions, which has several limitations. First, the final solution heavily
depends on heuristics or initial examples. Second, the number of LLM calls
might be unpredictable when the problem becomes complicated, such as TSP with
time window constraints. Third, Yang et al. [16], Liu et al. [17], and Liu et
al. [18] only investigate a single TSP problem. Ye et al. [19] only
investigate three variants of a single TSP problem, which is not enough to
demonstrate the ability of their methods on VRPs.

To overcome the limitations of these works, we construct a dataset with 21
variants of VRPs for single and multiple vehicles. We also propose a framework
that asks LLMs to generate solutions only with NL task descriptions.

Figure 3: Ten vehicle routing problems for a single robot.  Environment
information is context. Task is from users.

###  II-B LLMs for Code Generation and Code Repair

LLMs have shown remarkable progress in text-to-code generation without
specialized finetuning [21, 15, 22]. However, the generated code can sometimes
be flawed or incorrect [23]. There are various works that have been proposed
to fix this problem. We split them into two groups. The first group of works
aims to generate multiple solutions and then select the best solution, where
the criteria for determining ’the best’ differ. For example, Wang et al. [24]
define the best solution is the most consistent solution among all solutions;
Zhang et al. [25] ask LLMs to output a score for each solution and define the
best solution is the solution with the highest score. The second group of
works aims to generate one solution first and then ask LLMs to regenerate
solutions based on feedback (e.g., execution results). This process of
solution regeneration is known as self-reflection or self-refine. There are
various types of self-reflection given different feedback. Madaan et al. [26]
combine task-specific prompts with LLMs-generated outputs as feedback. Zhang
et al. [27] use the learned principles from previous incorrect solutions of
LLMs as feedback. Chen et al. [28] assume unit tests are given, combine code
explanation, execution results from the generated code, and unit tests as
feedback. However, high-quality feedback is often unavailable or heavily
relied on human effort [29].

In this work, we consider the fact that constraints in VRPs are naturally
described within the task—for example, a vehicle must visit all cities exactly
once. To amplify the benefits of utilizing these constraints, we first
instruct LLMs to extract constraints from NL task descriptions. Then, we ask
LLMs to generate unit tests given these constraints.

##  III Problem Formulation

As shown in Figure 2, given NL task descriptions of VRPs from users, we focus
on applying LLMs to perform text-to-Python code generation and then execute
the generated program to get the solution. However, to our knowledge, no VRP
dataset has been designed for LLMs. Specifically, no dataset contains various
VRPs with a list of (natural language task descriptions, solution) pairs.
Therefore, we construct a VRP dataset that contains ten variants of TSPs for a
single vehicle (Figure 3) and 11 variants of TSPs for multiple vehicles
(Figure 4).

The following two sections discuss our tasks for single- and multi-TSPs.
Throughout our discussion, we use the terms ’robot’ and ’vehicle’
interchangeably.

###  III-A TSP - Single Vehicle

For a single vehicle, we design the following ten tasks:

1) TSP, which aims to find the shortest tour while the robot needs to visit
each city exactly once. This is the default setting for the following tasks
unless otherwise specified. 2) Bottleneck TSP (BTSP), which aims to minimize
the maximum distance between any two consecutive cities in the tour. 3)
Generalized TSP (GTSP), also known as the traveling politician problem, aims
to visit each city group exactly once. 4) Sequential Ordering Problem (SOP) is
TSP with sequential ordering constraints. 5) K-TSP, which asks the vehicle to
visit selected cities only. 6) Multiple-Visit TSP (MV-TSP) allows the vehicle
to visit each city multiple times. 7) TSP with Time Windows (TSPTW) is TSP
with time window constraints. 8) Traveling Purchaser Problem (TPP) asks the
vehicle to visit a subset of cities to buy products while minimizing the total
cost (including travel cost and purchasing costs). 9) Multiple Depot TSP (MD-
TSP) without Capacity has multiple depots in the environment. 10) Multiple
Depot TSP (MD-TSP) with Capacity, different from 9), gives the vehicle travel
capacities, forcing it to visit a depot to recharge before it runs out of
energy.

Figure 4: Eleven vehicle routing problems for multiple robots.  Environment
information is context. Task is from users.

###  III-B TSP - Multiple Vehicle

For multiple vehicles, we design the following eleven tasks: 1) MTSP, where
all robots collaborate to visit all cities exactly once, aiming to minimize
the maximum tour cost among all robots. This is the default setting unless
specified otherwise. 2) Bottleneck MTSP (BMTSP) aims to minimize the maximum
distance between any two consecutive cities in the tour of all robots. 3)
Generalized MTSP (GMTSP) aims to minimize the maximum tour among all vehicles,
with each city group visited exactly once. 4) Sequential Ordering Problem
(MSOP) is MTSP with sequential ordering constraints. 5) K-MTSP, which asks
robots to visit selected cities only. 6) Multiple visit MTSP (MV-MTSP) allows
robots to visit each city multiple times. 7) MTSP with time windows (MTSPTW)
is MTSP with time window constraints. 8) Traveling Purchaser Problem (MTPP),
where all robots collaboratively visit a subset of cities to buy products,
aiming to minimize the total cost (including travel and purchasing costs). 9)
Multiple Depot MTSP (MD-MTSP) Without Capacity, multiple depots exist in the
environment. 10) Multiple depot MTSP (MD-MTSP) with Capacity, differing from
9), where all robots have travel capacities, necessitating visits to a depot
to recharge before they run out of energy. 11) Colored MTSP, where each city
has its color(s), and each robot has its own color. Each robot can only visit
cities that match its color.

##  IV Method

In this section, we start by searching the best basic prompt paradigms from
four basic prompt paradigms. We then propose a framework with self-reflection,
including self-debugging and self-verification, to correct the incorrect
programs or solutions generated at the first attempt.

###  IV-A Four Basic Prompt Paradigms

When addressing a vehicle routing problem described by natural language,
people usually translate the problem to mathematical formulations first and
then develop an algorithm to solve it. That is mainly because translating NL
task descriptions to mathematical formulations usually helps us to condense
the information and better understand the problem by filtering out unnecessary
information. This leads us to wonder: Would translating NL task descriptions
to intermediate descriptions, such as mathematical formulations, help LLMs
better solve VRPs? Moreover, when implementing algorithms to solve VRPs,
people usually use external solvers or libraries to assist in solving the
issues. Utilizing external solvers or libraries can simplify the
implementation of algorithms and enhance their efficiency compared to
developing algorithms from scratch and relying solely on Python’s standard
libraries. However, we observe that LLMs barely use external solvers or
libraries if we do not specify in the input prompt. This also leads us to
wonder: Would asking LLMs to solve VRPs with a specified external solver or
library help them better solve VRPs?

Based on the above two questions, we design four input prompt paradigms:

  * •

Generating Python code from NL task descriptions.

  * •

Translating NL task descriptions into mathematical formulations and then
generating Python code given mathematical formulations.

  * •

Generating Python code from NL task descriptions with the integration of a
specified external solver or library (e.g., Gurobi).

  * •

Translating NL task descriptions into mathematical formulations and then
generating Python code given mathematical formulations by integrating a
specified external solver or library (e.g., Gurobi).

###  IV-B Framework: Self-debugging and Self-verification

Inspired by the fact that high-quality external feedback can assist LLMs in
correcting their responses but high-quality external feedback (e.g., human
feedback) is often unavailable or heavily relied on human effort [29, 26, 28],
we propose a framework to generate unit tests based on constraints in the
problem descriptions. Unit tests are automated tests that ensure the program
meets its design and behaves as intended. Figure 2 shows an overview of our
framework.

In our framework, we ask LLMs to generate text-to-Python code based on NL task
descriptions from users. We then execute the Python code to get the execution
results. The execution results can be split into two situations. In the first
situation, the generated Python code might have bugs, and the execution may be
unsuccessful. In this circumstance, we concatenate problem description (task
description and environment information), Python code, and execution results
(i.e., traceback). We then send the concatenated prompt to LLMs and ask them
to regenerate the Python code by fixing bugs. In the second situation, the
generated Python code is executed successfully, and we get the solution. We
then generate unit tests to verify whether the solution is correct. We first
ask LLMs to extract constraints from the problem description. We then
concatenate problem descriptions, constraints, and execution results. After
that, we ask LLMs to generate unit tests given concatenated prompts. Finally,
we execute the unit tests to verify whether the solution is correct. If the
solution passes unit tests, we output the solution. Otherwise, we concatenate
the problem description, execution results, and verification analysis as
feedback and ask LLMs to regenerate Python code. We set the maximum number of
code reflections to 6.

##  V Experiments

We define three metrics to evaluate the performance. Feasibility: the ratio of
feasible solutions out of all solutions. Optimality: the ratio of optimal
solutions out of all solutions. Efficiency: the average of all normalized
solutions in $[0,1]$ by assigning $0$ to incorrect solutions and $1$ to
optimal solutions. We evaluate two LLMs, GPT-4 (gpt-4-0125-preview) and Gemini
1.0 Pro (version: gemini-1.0-pro-001).

In this work, our goal is to investigate the performance of LLMs in solving
vehicle routing problems. To this end, we seek to answer the following
questions.

###  V-A What is the best basic prompt paradigm among four basic prompt
paradigms?

Figure 5: The performance of GPT-4 in four basic prompt paradigms. We do not
use our proposed framework here. Each task is evaluated ten times separately.

In this section, we evaluate four basic prompt paradigms mentioned in Section
IV-A on GPT-4 and Gemini 1.0 Pro without using our proposed framework.

We first evaluate the performance of GPT-4. Figure 5 shows that, unlike humans
who solve VRPs, GPT-4 performs best with NL task descriptions without
specifying any external solvers or libraries. We analyze the generated code of
the prompt paradigms that use a specified external solver or library. These
prompt paradigms perform poorly because most of their generated programs have
bugs and are not executable. That shows GPT-4 faces challenges in generating
programs that rely on external solvers or libraries. We analyze the generated
programs and mathematical formulations of the prompt paradigm based on
mathematical formulations. We find that when GPT-4 translates NL to
mathematical formulations, there is a $\sim 50\%$ chance that it would discard
the position information of all cities, significantly decreasing the
performance of this prompt paradigm.

We then evaluate the performance of Gemini 1.0 Pro. We find more than $90\%$
of generated programs have bugs and are not executable, and the rest provide
incorrect solutions. Although there is a $1\%$ probability Gemini 1.0 Pro
providing the correct solution, we conclude that Gemini 1.0 Pro can not solve
VRPs at the initial attempt.

Figure 6: The performance of GPT-4 in four basic prompt paradigms. We use our
proposed framework and set the maximum reflection number to six. Each task is
evaluated ten times separately.

###  V-B What is the performance of LLMs with our proposed framework?

In this section, we evaluate the performance of LLMs with our proposed
framework across four basic prompt paradigms mentioned in Section IV-A. We
conduct four studies.

First, we evaluate the performance of GPT-4 in four paradigms with our
proposed framework. Figure 6 shows the prompt paradigm that directly uses NL
still performs the best.

Second, we evaluate the efficacy of our proposed framework on GPT-4. Figure 7
shows that our proposed framework can steadily improve the performance of
GPT-4 with the increase of reflection times. Table I shows that with our
proposed framework when the maximum reflection number approaches 3, our
framework can improve $10\%$ feasibility, $4\%$ optimality, and $9\%$
efficiency compared with without our proposed framework (maximum reflection
number is 0). When the maximum reflection number approaches 6, our framework
can further improve $6\%$ feasibility, $3\%$ optimality, and $6\%$ efficiency.
Besides, we have three interesting findings. First, the performance of GPT-4
in the environment with one vehicle and five cities is worse than in the
environment with one vehicle and ten cities. We hypothesize that GPT-4 is
inconsistent in solving problems. While a performance discrepancy is evident
between scenarios of one vehicle with ten cities and four vehicles with ten
cities, such discrepancies become inconsistent in similar difficulty levels,
such as between environments of one vehicle with five and ten cities. Second,
task 3-GTSP shows that reflection negatively impacts performance in the one-
vehicle and ten-city scenarios, as indicated by the optimality chart. We
hypothesize that GPT-4-generated unit tests sometimes misclassify correct
solutions as incorrect. Third, as the number of vehicles increases to four,
the performance of GPT-4 significantly decreases. Although the performance of
GPT-4 increases with our proposed framework, the performance can not match the
scenario with one vehicle and ten cities. This illustrates that GPT-4
struggles to effectively manage scenarios involving four vehicles and ten
cities. As tasks grow increasingly complex, self-reflection reveals its
limitations, failing to match the performance of LLMs in situations where the
difficulty levels of the tasks vary.

Figure 7: The relationship between performance and the number of reflections.
In the legend, $0$ means no reflection, $3$ means the maximum reflection
number is 3, and $6$ means the maximum reflection number is 6. Each task is
evaluated ten times separately.

Third, we evaluate the performance of Gemini 1.0 Pro in four paradigms with
our proposed framework and set the maximum reflection number to six. However,
the performance of Gemini 1.0 Pro does not improve in any prompt paradigms.
That shows Gemini 1.0 Pro cannot fix bugs in its previously generated programs
or correct its previous incorrect solutions with our proposed framework.

TABLE I: Relationship Between Performance and Reflection Number

| Feasibility $\uparrow$ | Optimality $\uparrow$ | Efficiency $\uparrow$
---|---|---|---
| 1-5 | 1-10 | 4-10 | all | 1-5 | 1-10 | 4-10 | all | 1-5 | 1-10 | 4-10 | all
0 | 0.69 | 0.75 | 0.23 | 0.56 | 0.48 | 0.65 | 0.08 | 0.40 | 0.65 | 0.73 | 0.20 | 0.53
3 | 0.80 | 0.83 | 0.35 | 0.66 | 0.50 | 0.68 | 0.13 | 0.44 | 0.74 | 0.80 | 0.32 | 0.62
6 | 0.83 | 0.90 | 0.44 | 0.72 | 0.51 | 0.75 | 0.16 | 0.47 | 0.77 | 0.87 | 0.39 | 0.68

###  V-C Can adding more details to task descriptions increase the performance
of GPT-4?

When we ask LLMs to solve problems, we usually struggle with how to design
thorough task descriptions. One question is: Would the performance of LLMs
drop if task descriptions were without some details that further explain the
task, despite the core meaning of tasks being intact from a human perspective?
To explore this, we simplify our task descriptions by removing details that
further explain tasks, ensuring the core meaning remains intact. We apply this
approach to all tasks in our dataset and create a new dataset. We then
evaluate the performance of GPT-4 on this new dataset, using our framework and
a prompt paradigm that instructs GPT-4 to generate code based on NL
descriptions. Table II shows that after removing details that further explain
tasks, the performance of GPT-4 drops $4\%$ in feasibility, $4\%$ in
optimality, and $5\%$ in efficiency. The results demonstrate the importance of
adding details to further explain tasks. However, although details that
further explain tasks can increase the performance of GPT-4 in general, Figure
8 shows that not all tasks can benefit from these details, and these details
might harm the performance of GPT-4. We hypothesize that adding details to
further explain tasks might make task descriptions redundant, potentially
hindering the ability of GPT-4 to understand tasks.

Figure 8: Comparing the performance of GPT-4 in task descriptions with and
without details to explain tasks further. Each task is evaluated ten times
separately. TABLE II: W/O Details vs W/ Details

| Feasibility $\uparrow$ | Optimality $\uparrow$ | Efficiency $\uparrow$
---|---|---|---
| 1-5 | 1-10 | 4-10 | all | 1-5 | 1-10 | 4-10 | all | 1-5 | 1-10 | 4-10 | all
W/O Details | 0.80 | 0.82 | 0.42 | 0.68 | 0.59 | 0.59 | 0.12 | 0.43 | 0.76 | 0.78 | 0.35 | 0.63
W/ Details | 0.83 | 0.90 | 0.44 | 0.72 | 0.51 | 0.75 | 0.16 | 0.47 | 0.77 | 0.87 | 0.39 | 0.68

###  V-D Identifying what details can benefit the performance of LLMs is
challenging. Is there a method that can help us?

In the above question, we find that details that further explain tasks may or
may not decrease the performance of LLMs. However, identifying if adding more
information to the initial task descriptions can benefit the performance of
LLMs is challenging.

In this section, we propose a mechanism that allows LLMs to identify if task
descriptions are thorough, propose questions to users if they think task
descriptions are not thorough, and rephrase task descriptions by absorbing
users’ answers into initial task descriptions. We limit the number of LLMs
requesting assistance to three. Figure 9 shows an overview of our mechanism.

Figure 9: Mechanism overview: ask users for help

We evaluate the performance of this mechanism by selecting three typical tasks
(the performance gap between task descriptions with and without details is
biggest across all tasks in Figure 8) from three different environment
settings on GPT-4. We select the SOP task from the environment with one
vehicle and five cities. We select the TSPTW task from the environment with
one vehicle and ten cities. We select the CMTSP task from the environment with
four vehicles and ten cities. Table III shows that our mechanism can increase
the performance of GPT-4 in general. However, the performance increase varies
for different tasks. For the SOP task in the environment with one vehicle and
ten cities, the ratio of asking for help from GPT-4 is only three out of ten.
This indicates that GPT-4 struggles to efficiently identify missing details,
consequently limiting its performance. For the TSPTW task in the environment
with one vehicle and ten cities, the ratio of asking for help from GPT-4 is
eight out of ten times, showing that GPT-4 can identify the task as not
thorough and ask questions to make task descriptions thorough. For the CMTSP
task in the environment with four vehicles and ten cities, the ratio of asking
for help from GPT-4 is nine out of ten times. However, the performance of
GPT-4 does not increase. We dive into the questions that GPT-4 asks and find
that GPT-4 keeps asking unnecessary questions (e.g., How many cities are
there?). This finding shows that although GPT-4 considers that task
descriptions are not thorough, it does not correctly identify the needed
information.

TABLE III: W/O Mechanism vs W/ Mechanism

| SOP: 1 vehicle, 5 cities | TSPTW: 1 vehicle, 10 cities | CMTSP: 4 vehicles, 10 cities
---|---|---|---
| Feasible $\uparrow$ | Optimal $\uparrow$ | Efficient $\uparrow$ | Ratio of ask $\uparrow$ | Feasible $\uparrow$ | Optimal $\uparrow$ | Efficient $\uparrow$ | Ratio of ask $\uparrow$ | Feasible $\uparrow$ | Optimal $\uparrow$ | Efficient $\uparrow$ | Ratio of ask $\uparrow$
W/O Mechanism | 0.50 | 0.00 | 0.40 | - | 0.50 | 0.50 | 0.50 | - | 0.40 | 0.00 | 0.29 | -
W/ Mechanism | 0.50 | 0.30 | 0.46 | 3/10 | 0.90 | 0.90 | 0.90 | 8/10 | 0.40 | 0.00 | 0.31 | 9/10

###  V-E What is the performance of Claude 3 in VRPs?

The recent promising performance of Claude 3 leads us to wonder: What is the
performance of Claude 3 in VRPs? Will the performance of Claude 3 outperform
GPT-4? We integrate Claude 3 with our proposed framework and evaluate the
performance of Claude 3 by directly feeding NL into Claude 3. We select
claude-3-opus-20240229, the most powerful model of Claude. Table IV shows that
although Claude 3 Opus can solve VRPs, its performance is worse than GPT-4.

TABLE IV: Claude 3 vs GPT-4

| Feasibility $\uparrow$ | Optimality $\uparrow$ | Efficiency $\uparrow$
---|---|---|---
| 1-5 | 1-10 | 4-10 | all | 1-5 | 1-10 | 4-10 | all | 1-5 | 1-10 | 4-10 | all
Claude 3 | 0.47 | 0.58 | 0.30 | 0.45 | 0.29 | 0.11 | 0.07 | 0.16 | 0.43 | 0.49 | 0.29 | 0.40
GPT-4 | 0.83 | 0.90 | 0.44 | 0.72 | 0.51 | 0.75 | 0.16 | 0.47 | 0.77 | 0.87 | 0.39 | 0.68

##  VI Conclusion and Future Work

In this work, we investigate the ability of Large Language Models to solve
Vehicle Routing Problems. We start by constructing a dataset with 21 VRPs for
single and multiple vehicles. We search four basic prompt paradigms for LLMs
and find that directly feeding natural language into LLMs provides the best
performance. Next, we propose a framework with self-reflection to improve the
performance of LLMs in solving VRPs. Additionally, we study the importance of
details in explaining tasks and propose a mechanism to help LLMs get thorough
task descriptions. We finally evaluate the performance of Claude 3 Opus on
VRPs and find the performance of Claude 3 Opus is better than Gemini 1.0 Pro
but worse than GPT-4.

## References

  * [1] H. Wu, F. Tao, and B. Yang, “Optimization of vehicle routing for waste collection and transportation,” _International Journal of Environmental Research and Public Health_ , vol. 17, no. 14, p. 4963, 2020.
  * [2] D. Prajapati, F. T. Chan, Y. Daultani, and S. Pratap, “Sustainable vehicle routing of agro-food grains in the e-commerce industry,” _International Journal of Production Research_ , vol. 60, no. 24, pp. 7319–7344, 2022.
  * [3] A. Botros, B. Gilhuly, N. Wilde, A. Sadeghi, J. Alonso-Mora, and S. L. Smith, “Optimizing task waiting times in dynamic vehicle routing,” _IEEE Robotics and Automation Letters_ , 2023.
  * [4] L. Wenyi, “Research on vehicle routing problem and application scenarios,” in _Handbook of Mobility Data Mining_. Elsevier, 2023, pp. 63–88.
  * [5] J. Mańdziuk and A. Żychowski, “A memetic approach to vehicle routing problem with dynamic requests,” _Applied Soft Computing_ , vol. 48, pp. 522–534, 2016.
  * [6] K. K. Ng, C. K. Lee, S. Zhang, K. Wu, and W. Ho, “A multiple colonies artificial bee colony algorithm for a capacitated vehicle routing problem and re-routing strategies under time-dependent traffic congestion,” _Computers & Industrial Engineering_, vol. 109, pp. 151–168, 2017.
  * [7] D. Zhang, S. Cai, F. Ye, Y.-W. Si, and T. T. Nguyen, “A hybrid algorithm for a vehicle routing problem with realistic constraints,” _Information Sciences_ , vol. 394, pp. 167–182, 2017.
  * [8] A. Rashid, S. Sharma, C. M. Kim, J. Kerr, L. Y. Chen, A. Kanazawa, and K. Goldberg, “Language embedded radiance fields for zero-shot task-oriented grasping,” in _7th Annual Conference on Robot Learning_ , 2023.
  * [9] A. Z. Ren, A. Dixit, A. Bodrova, S. Singh, S. Tu, N. Brown, P. Xu, L. Takayama, F. Xia, J. Varley _et al._ , “Robots that ask for help: Uncertainty alignment for large language model planners,” _arXiv preprint arXiv:2307.01928_ , 2023.
  * [10] J. Zhang, J. Zhang, K. Pertsch, Z. Liu, X. Ren, M. Chang, S.-H. Sun, and J. J. Lim, “Bootstrap your own skills: Learning to solve new tasks with large language model guidance,” _arXiv preprint arXiv:2310.10021_ , 2023.
  * [11] K. Rana, J. Haviland, S. Garg, J. Abou-Chakra, I. Reid, and N. Suenderhauf, “Sayplan: Grounding large language models using 3d scene graphs for scalable task planning,” _arXiv preprint arXiv:2307.06135_ , 2023.
  * [12] D. Shah, M. R. Equi, B. Osiński, F. Xia, B. Ichter, and S. Levine, “Navigation with large language models: Semantic guesswork as a heuristic for planning,” in _Conference on Robot Learning_. PMLR, 2023, pp. 2683–2699.
  * [13] Y. Tang, W. Yu, J. Tan, H. Zen, A. Faust, and T. Harada, “Saytap: Language to quadrupedal locomotion,” _arXiv preprint arXiv:2306.07580_ , 2023.
  * [14] Y. Chen, J. Arkin, Y. Zhang, N. Roy, and C. Fan, “Autotamp: Autoregressive task and motion planning with llms as translators and checkers,” _arXiv preprint arXiv:2306.06531_ , 2023.
  * [15] E. Nijkamp, B. Pang, H. Hayashi, L. Tu, H. Wang, Y. Zhou, S. Savarese, and C. Xiong, “Codegen: An open large language model for code with multi-turn program synthesis,” _arXiv preprint arXiv:2203.13474_ , 2022.
  * [16] C. Yang, X. Wang, Y. Lu, H. Liu, Q. V. Le, D. Zhou, and X. Chen, “Large language models as optimizers,” _arXiv preprint arXiv:2309.03409_ , 2023.
  * [17] S. Liu, C. Chen, X. Qu, K. Tang, and Y.-S. Ong, “Large language models as evolutionary optimizers,” _arXiv preprint arXiv:2310.19046_ , 2023.
  * [18] F. Liu, X. Tong, M. Yuan, and Q. Zhang, “Algorithm evolution using large language model,” _arXiv preprint arXiv:2311.15249_ , 2023.
  * [19] H. Ye, J. Wang, Z. Cao, and G. Song, “Reevo: Large language models as hyper-heuristics with reflective evolution,” _arXiv preprint arXiv:2402.01145_ , 2024.
  * [20] Y. Huang, W. Zhang, L. Feng, X. Wu, and K. C. Tan, “How multimodal integration boost the performance of llm for optimization: Case study on capacitated vehicle routing problems,” _arXiv preprint arXiv:2403.01757_ , 2024.
  * [21] J. Li, B. Hui, R. Cheng, B. Qin, C. Ma, N. Huo, F. Huang, W. Du, L. Si, and Y. Li, “Graphix-t5: Mixing pre-trained transformers with graph-aware layers for text-to-sql parsing,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 37, no. 11, 2023, pp. 13 076–13 084.
  * [22] Q. Zheng, X. Xia, X. Zou, Y. Dong, S. Wang, Y. Xue, Z. Wang, L. Shen, A. Wang, Y. Li _et al._ , “Codegeex: A pre-trained model for code generation with multilingual evaluations on humaneval-x,” _arXiv preprint arXiv:2303.17568_ , 2023.
  * [23] L. Pan, M. Saxon, W. Xu, D. Nathani, X. Wang, and W. Y. Wang, “Automatically correcting large language models: Surveying the landscape of diverse self-correction strategies,” _arXiv preprint arXiv:2308.03188_ , 2023.
  * [24] X. Wang, J. Wei, D. Schuurmans, Q. Le, E. Chi, S. Narang, A. Chowdhery, and D. Zhou, “Self-consistency improves chain of thought reasoning in language models,” _arXiv preprint arXiv:2203.11171_ , 2022.
  * [25] T. Zhang, T. Yu, T. Hashimoto, M. Lewis, W.-t. Yih, D. Fried, and S. Wang, “Coder reviewer reranking for code generation,” in _International Conference on Machine Learning_. PMLR, 2023, pp. 41 832–41 846.
  * [26] A. Madaan, N. Tandon, P. Gupta, S. Hallinan, L. Gao, S. Wiegreffe, U. Alon, N. Dziri, S. Prabhumoye, Y. Yang _et al._ , “Self-refine: Iterative refinement with self-feedback,” _Advances in Neural Information Processing Systems_ , vol. 36, 2024.
  * [27] T. Zhang, A. Madaan, L. Gao, S. Zheng, S. Mishra, Y. Yang, N. Tandon, and U. Alon, “In-context principle learning from mistakes,” _arXiv preprint arXiv:2402.05403_ , 2024.
  * [28] X. Chen, M. Lin, N. Schärli, and D. Zhou, “Teaching large language models to self-debug,” _arXiv preprint arXiv:2304.05128_ , 2023.
  * [29] J. Huang, X. Chen, S. Mishra, H. S. Zheng, A. W. Yu, X. Song, and D. Zhou, “Large language models cannot self-correct reasoning yet,” _arXiv preprint arXiv:2310.01798_ , 2023.