\section{Analysis of sensitivity to inputs}
\label{sec:sens}

<<echo=false,results=hide>>=
library(tgp)
seed <- 0; set.seed(seed)
@ 


%% Bobby says to Taddy:  Need to check that LHS stuff does the right
%% thing for the inputs that are being treated as categorical.  Not
%% sure what this means because exactly since "treating as categorical"
%% is means "no-GP" or "tree-only" on those inputs.  One way might be
%% to allow a scale=0 argument at which point the resulting LHS samples
%% snap to {0,1}  The mode argument could determine the p parameter
%% in the Bernoulli(p) distribution that would result.
%% Taddy replies: done and done!  At least, the following seems fine.
%% re-fit the ozone data with temperature as a threshold effect:
%% library(tgp)
%% X <- airquality[,2:4]
%% X$Temp[X$Temp >74] <- 1
%% X$Temp[X$Temp >1] <- 0
%% Z <- airquality$Ozone
%% rect <- t(apply(X, 2, range, na.rm=TRUE))
%% mode <- apply(X , 2, mean, na.rm=TRUE)
%% shape <- rep(2,3)
%% mode[3] <- .5
%% shape[3] <- 0
%% s.air <- sens(X=X, Z=Z, nn.lhs=300, rect=rect, 
%% shape=shape, mode=mode, basemax=2)
%% plot(s.air, layout="sens")
%% plot(s.air, layout="sens", maineff=c(1,2,3))
%%
%% This example has been added to the vignette.
%% This won't work for your 'binary' representation of
%% categorical variables with more than two levels, however.
%% It will provide an answer, but perhaps not what you'd want
%% since it will look at the sensitivity to each binary element,
%% rather than at sensitivity to the unique sequences. 
%% This is fixable, but it would take some real effort and time (I'll make
%% sure to get it done for a future iteration).
%% The code for plotting is a bit messy, but it should be robust enough.
%% I've added a note in the description of 'shape' in sens.Rd outlining
%% these changes.
%% Changes in the c-code only occur in tgp.cc, where the
%% main effects calculation now takes basemax into consideration.

Methods for understanding how inputs, or explanatory variables,
contribute to the outputs, or response, of simple statistical models
are by now classic in the literature and frequently used in
practical application.  For example, in linear regression one can
perform $F$--tests to ascertain the relevance of a predictor, or inspect
the leverage of a particular input setting, or use Cooks' distance, to
name a few.  Unfortunately, such convenient statistics/methods are not
available for more complicated models, such as those in the {\tt tgp}
family of nonparametric models.  A more advanced tool is needed.

Sensitivity Analysis (SA) is a resolving of the sources of output
variability by apportioning elements of this variation to different
sets of input variables.  It is applicable in wide generality.  The
edited volume by Saltelli et al.~\cite{SaltChanScot2000} provides an
overview of the field.  Valuable recent work on smoothing methods is
found in \cite{StorHelt2008,VeigWahlGamb2009}, and Storlie, et
al.~\cite{StorSwilHeltSall2009}, provide a nice overview of
nonparametric regression methods for inference about sensitivity. The
analysis of response variability is useful in a variety of different
settings.  For example, when there is a large number of input
variables over which an objective function is to be optimized,
typically only a small subset will be influential within the confines
of their uncertainty distribution.  SA can be used to reduce the input
space of such optimizations \cite{TaddLeeGrayGrif2009}.  Other authors
have used SA to assess the risk associated with dynamic factors
affecting the storage of nuclear waste \cite{HommSalt1996}, and to
investigate the uncertainty characteristics of a remote sensing model
for the reflection of light by surface vegetation
\cite{MorrKottTaddFurfGana2008}.  The {\tt sens} function adds to {\tt
  tgp} a suite of tools for global sensitivity analysis, and enables
``out-of-the-box'' estimation of valuable sensitivity indices for any
regression relationship that may be modeled by a member of the {\tt
  tgp} family.

The type of sensitivity analysis provided by {\tt tgp} falls within
the paradigm of global sensitivity analysis, wherein the variability
of the response is investigated with respect to a probability
distribution over the entire input space.  The recent book by Saltelli
et al. \cite{SaltEtAl2008} serves as a primer on this field.  Global
SA is inherently a problem of statistical inference, as evidenced by
the interpolation and estimation required in a study of the full range
of inputs.  This is in contrast with the analytical nature of local
SA, which involves derivative--based investigation of the stability of
the response over a small region of inputs.  We will ignore local SA
for the remainder of this document.

The sensitivity of a response $z$ to a changing input $\mb{x}$ is
always considered in relation to a specified {\it uncertainty
  distribution}, defined by the density $u(\mb{x})$, and the
appropriate marginal densities $u_i(x_i)$.  What is represented by the
uncertainty distribution changes depending upon the context.  The
canonical setup has that $z$ is the response from a complicated
physics or engineering simulation model, with tuning parameters
$\mb{x}$, that is used to predict physical phenomena. In this
situation, $u(\mb{x})$ represents the experimentalist's uncertainty
about real--world values of $\mb{x}$.  In optimization problems, the
uncertainty distribution can be used to express prior information from
experimentalists or modelers on where to look for solutions.  Finally,
in the case of observational systems (such as air-quality or smog
levels), $u(\mb{x})$ may be an estimate of the density governing the
natural occurrence of the $\mb{x}$ factors (e.g., air-pressure,
temperature, wind, and cloud cover).  In this setup, SA attempts to
resolve the natural variability of $z$.

The most common notion of sensitivity is tied to the relationship
between conditional and marginal variance for $z$.  Specifically,
variance--based methods decompose the variance of the objective
function, with respect to the uncertainty distribution on the inputs,
into variances of conditional expectations.  These are a natural
measure of the output association with specific sets of variables and
provide a basis upon which the importance of individual inputs may be
judged.  The other common component of global SA is an accounting of
the main effects for each input variable, $\mathbb{E}_{u_j}[z|x_j]$,
which can be obtained as a by-product of the variance analysis.  

Our variance--based approach to SA is a version of the method of
Sobol', wherein a deterministic objective function is decomposed into
summands of functions on lower dimensional subsets of the input
space. Consider the function decomposition $ f(x_1, \ldots ,x_d) = f_0
+ \sum_{j=1}^df_j(x_j) +\sum_{1 \leq i < j \leq d} f_{ij}(x_j,x_i) +
\ldots + f_{1,\ldots,d}(x_1, \ldots ,x_d).  $ When the response $f$ is
modeled as a stochastic process $z$ conditional on inputs $\mb{x}$, we
can develop a similar decomposition into the response distributions
which arise when $z$ has been marginalized over one subset of
covariates and the complement of this subset is allowed to vary
according to a marginalized uncertainty distribution.  In particular,
we can obtain the marginal conditional expectation
$\mbb{E}[z|\mb{x}_J=\{x_j:j\in J\}]$ $=$ $\int_{\mathbb{R}^{d-d_J}}
\mbb{E}[z|\mb{x}]u(\mb{x}) d\mb{x}_{-J}$, where $J=\{j_1, \ldots,
j_{d_J}\}$ indicates a subset of input variables, $\mb{x}_{-j}
=\{x_j:j\notin J\}$, and the marginal uncertainty density is given by
$u_J(\mb{x}_J) = \int_{\mathbb{R}^{d-d_J}} u(\mb{x}) d\{x_i:i \notin J
\}$.  SA concerns the variability of $\mbb{E}[z|\mb{x}_J]$ with
respect to changes in $\mb{x}_J$ according to $u_J(\mb{x}_J)$ and, if
$u$ is such that the inputs are uncorrelated, the variance
decomposition is available as
\begin{equation} \label{eqn:var_decomp}
  \mr{var}(\mbb{E}[z|\mb{x}]) = \sum_{j=1}^dV_j
+ \sum_{1 \leq i < j \leq d} V_{ij} + \ldots + V_{1,\ldots,d},
\end{equation}
where $V_j = \mr{var}(\mbb{E}[z|x_j])$,
$V_{ij}=\mr{var}(\mbb{E}[z|x_i, x_j]) - V_i - V_j$, and so on.   Clearly, when the
inputs are correlated this identity no longer holds (although a
``less-than-or-equal-to'' inequality is always true).  But it is
useful to retain an intuitive interpretation of the $V_J$'s as a
portion of the overall marginal variance.

Our global SA will focus on the related sensitivity indices $S_J =
V_J/\mr{var}(z)$ which, as can be seen in the above equation, will
sum to one over all possible $J$ and are bounded to $[0,1]$.  These
$S_J$'s provide a natural measure of the {\it importance} of a set $J$
of inputs and serve as the basis for an elegant analysis of
sensitivity.  The {\tt sens} function allows for easy calculation of
two very important sensitivity indices associated with each input:
the 1$^{\rm st}$ order for the $j$th input variable,
\begin{equation}
S_j = \frac{\mr{var}\left(\mbb{E}\left[z|x_j\right]\right)}{\mr{var}(z)},
\label{eq:S}
\end{equation}
and the total sensitivity for input $j$, 
\begin{equation}
T_j = \label{eq:T} \frac{\mbb{E}\left[\mr{var}\left(z|\mb{x}_{-j}\right)\right]}{\mr{var}(z)}.
\end{equation}
The 1$^{\rm st}$ order indices measure the portion of variability that is due
to variation in the main effects for each input variable, while the
total effect indices measure the portion of variability that is due to
total variation in each input.  From the identity
$\mbb{E}\left[\mr{var}\left(z|\mb{x}_{-j}\right)\right] = \mr{var}(z)
- \mr{var}\left(\mbb{E}\left[z|\mb{x}_{-j}\right]\right)$, it can be
seen that $T_j$ measures the {\it residual} variability remaining
after variability in all other inputs has been apportioned and that,
for a deterministic response and uncorrelated input variables, $T_j =
\sum_{J:j \in J} S_J$.  This implies that the difference between $T_j$
and $S_j$ provides a measure of the variability in $z$ due to
interaction between input $j$ and the other input variables.  A large
difference may lead the investigator to consider other sensitivity
indices to determine where this interaction is most influential, and
this is often a key aspect of the dimension--reduction that SA provides
for optimization problems.

\subsection{Monte Carlo integration for sensitivity indices}

Due to the many integrals involved, estimation of the sensitivity
indices is not straightforward.  The influential paper by Oakley \&
O'Hagan \cite{OaklOhag2004} describes an empirical Bayes estimation
procedure for the sensitivity indices, however some variability in the
indices is lost due to plug-in estimation of GP model parameters and,
more worryingly, the variance ratios are only possible in the form of
a ratio of expected values.  Marrel, et
al.~\cite{MarrIoosLaurRous2009}, provide a more complete analysis of
the GP approach to this problem, but their methods remain restricted
to estimation of the first order Sobol indices. Likelihood based
approaches have also been proposed
\cite{WelcBuckSackWynnMitcMorr1992,MorrKottTaddFurfGana2008}.  The
technique implemented in {\tt tgp} is, in contrast, fully Bayesian and
provides a complete accounting of the uncertainty involved.  Briefly,
at each iteration of an MCMC chain sampling from the treed GP
posterior, output is predicted over a large (carefully chosen) set of
input locations.  Conditional on this predicted output, the
sensitivity indices can be calculated via Monte Carlo integration.  By
conditioning on the predicted response (and working as though it were
the observed response), we obtain a posterior sample of the indices,
incorporating variability from both the integral estimation and
uncertainty about the function output.  In particular, the {\tt sens}
function includes a {\tt model} argument which allows for SA based on
any of the prediction models (the {\tt b*} functions) in {\tt tgp}.

Our Monte Carlo integration is based upon Saltelli's \cite{Salt2002}
efficient Latin hypercube sampling (LHS) scheme for estimation of both
1$^{\rm st}$ order and total effect indices. We note that the estimation is
only valid for uncorrelated inputs, such that $u(\mb{x}) =
\prod_{j=1}^d u_j(x_j)$.  The {\tt sens} function only allows for
uncertainty distributions of this type (in fact, the marginal
distributions also need to be bounded), but this is a feature of nearly
every ``out-of-the-box'' approach to SA.  Studies which concern
correlated inputs will inevitably require modeling for this
correlation, whereas most regression models (including those in {\tt
  tgp}) condition on the inputs and ignore the joint density for
$\mb{x}$.  Refer to the work of Saltelli \& Tarantola
\cite{SaltTara2002} for an example of SA with correlated inputs.


We now briefly describe the integration scheme.  The 2nd moment is a
useful intermediate quantity in variance estimation, and we define
\[
D_J = \mbb{E}\left[\mbb{E}^2\left[z|\mb{x}_{J}\right]\right] =
\int_{\mbb{R}^{d_J}} \mbb{E}^2\left[z|
  {\mb{x}_J}\right]u_J(\mb{x}_J)d(\mb{x}_J).
\]
Making use of an auxiliary variable,
\begin{eqnarray*}
  D_J &=&  \int_{\mbb{R}^{d_J}} \left[\int_{\mbb{R}^{d_{-J}}} 
    \!\!\!\mbb{E}\left[ z | \mb{x}_J, \mb{x}_{-J} \right]u_{-J}(\mb{x}_{-J})d\mb{x}_{-J}
    \int_{\mbb{R}^{d_{-J}}} \!\!\!\mbb{E}\left[ z | \mb{x}_J, \mb{x}'_{-J} \right]
    u_{-J}(\mb{x}'_{-J})d\mb{x}'_{-J}\right]u_J(\mb{x}_J)\mb{x}_{J}\\
  &=& \int_{\mbb{R}^{d + d_{-J}}} 
  \!\!\mbb{E}\left[ z | \mb{x}_J, \mb{x}_{-J} \right]\mbb{E}\left[ z | \mb{x}_J, \mb{x}'_{-J} \right] 
  u_{-J}(\mb{x}_{-J})u_{-J}(\mb{x}'_{-J})u_{J}(\mb{x}_{J})d\mb{x}d\mb{x}'_{J}. 
\end{eqnarray*}
Thus, in the case of independent inputs, 
\[ 
D_J =
\int_{\mbb{R}^{d+d_{-J}}} \mbb{E}\left[ z |\mb{x} \right]\mbb{E}\left[
  z | \mb{x}_J, \mb{x}'_{-J} \right] u_{-J}(\mb{x}'_{-J})u({\bf
x})d\mb{x}'_{-J}d\mb{x}.
\]  
Note that at this point, if the inputs had
been correlated, the integral would have been instead with respect to
the joint density $u(\mb{x})u(\mb{x}_{-J}' | \mb{x}_J)$, leading to a more
difficult integral estimation problem.

Recognizing that $S_j = (D_j-\mbb{E}^2[z])/\mr{var}(z)$ and $T_j = 1-
\left( \left(D_{-j} - \mbb{E}^2[z]\right)\right)/\mr{var}(z)$, we need
estimates of $\mr{var}(z)$, $\mbb{E}^2[z]$, and $\{ (D_j, D_{-j}) :
j=1,\ldots,d \}$ to calculate the sensitivity indices.  Given a LHS
$M$ proportional to $u(\mb{x})$,
\begin{equation*}
M = \left[ 
\begin{array}{c}
s_{1_1} ~ \cdots ~ s_{1_d}\\ 
\vdots  \\
s_{m_1} ~ \cdots ~ s_{m_d}\\
\end{array}
\right],
\end{equation*}
it is possible to estimate $\widehat{\mbb{E}[z]} = \frac{1}{m}
\sum_{k=1}^m\mbb{E}[z|{\bf s}_k]$ and $\widehat{\mr{var}[z]} =
\frac{1}{m} \mbb{E}^T[z|M]\mbb{E}[z|M] - \widehat{\mbb{E}[z]}\widehat{\mbb{E}[z]}$,
where the convenient notation $\mbb{E}[z|M]$ is taken to mean
$\left[\mbb{E}[z|\mb{s}_1] \cdots \mbb{E}[z|\mb{s}_m]\right]^T$.  All
that remains is to estimate the $D$'s.  Define a second LHS $M'$
proportional to $u$ of the same size as $M$ and say that $N_J$ is $M'$
with the $J$ columns replaced by the corresponding columns of $M$.
Hence,
\begin{equation*}
N_j = \left[ 
\begin{array}{c}
s'_{1_1} \cdots s_{1_j} \cdots s'_{1_d}\\ 
\vdots  \\
s'_{m_1} \cdots s_{m_j} \cdots s'_{m_d}
\end{array}\right]~~~\mr{and}~~~
N_{-j} = \left[ 
\begin{array}{c}
s_{1_1} \cdots s'_{1_j} \cdots s_{1_d}\\ 
\vdots  \\
s_{m_1} \cdots s'_{m_j} \cdots s_{m_d} 
\end{array}\right].
\end{equation*}
The estimates are then $\hat D_j =
\mbb{E}^T[z|M]\mbb{E}[z|N_{j}]/(m-1)$ and $\hat D_{-j}$ $=$
$\mbb{E}^T[z|M']\mbb{E}[z|N_{j}]/(m-1)$ $\approx$ $
\mbb{E}^T[z|M]\mbb{E}[z|N_{-j}]/(m-1)$. Along with the variance and
expectation estimates, these can be plugged into equations for $S_j$
and $T_j$ in (\ref{eq:S}--\ref{eq:T}) to obtain $\hat S_j$ and $\hat
T_j$.  Note that Saltelli recommends the use of the alternative
estimate $\widehat{\mbb{E}^2[z]} =
\frac{1}{n-1}\mbb{E}^T[z|M]\mbb{E}[z|M']$ in calculating 1$^{\rm st}$
order indices, as this brings the index closer to zero for
non-influential variables. However, it has been our experience that
these biased estimates can be unstable, and so {\tt tgp} uses the
standard $\widehat{\mbb{E}^2[z]} =
\widehat{\mbb{E}[z]}\widehat{\mbb{E}[z]}$ throughout.  As a final
point, we note that identical MCMC sampling-based integration schemes
can be used to estimate other Sobol indices (e.g., second order, etc)
for particular combinations of inputs, but that this would require
customization of the {\tt tgp} software.

The set of input locations which need to be evaluated for each
calculation of the indices is $\{ M, M', N_1,\ldots,N_d \}$, and if
$m$ is the sample size for the Monte Carlo estimate this scheme
requires $m(d+2)$ function evaluations. Hence, at each MCMC
iteration of the model fitting, the $m(d+2)$ locations are drawn
randomly according the LHS scheme, creating a random prediction
matrix, {\tt XX}.  By allowing random draws of the
input locations, the Monte Carlo error of the integral estimates will
be included in the posterior variability of the indices and the
posterior moments will not be dependent upon any single estimation
input set.  Using predicted output over this input set, a single
realization of the sensitivity indices is calculated through
Saltelli's scheme.  At the conclusion of the MCMC, we have a
representative sample from the posterior for ${\bf S}$ and ${\bf T}$.
The averages for these samples are unbiased estimates of the
posterior mean, and the variability of the sample is representative of
the complete uncertainty about model sensitivity.

Since a subset of the predictive locations ($M$ and $M'$) are actually
a LHS proportional to the uncertainty distribution, we can also
estimate the main effects at little extra computational cost. At each
MCMC iteration, a one--dimensional nonparametric regression is fit
through the scatterplot of $[s_{1_j}, \ldots, s_{m_j},s'_{1_j},
\ldots, s'_{m_j}]$ vs. $[\mbb{E}[z|M],\mbb{E}[z|M']]$ for each of the
$j=1,\ldots,d$ input variables.  The resultant regression estimate
provides a realization of $\mbb{E}[z|x_j]$ over a grid of $x_j$
values, and therefore a posterior draw of the main effect curve.
Thus, at the end of the MCMC, we have not only unbiased estimates of
the main effects through posterior expectation, but also a full
accounting of our uncertainty about the main effect curve.  This
technique is not very sensitive to the method of non-parametric
regression, since $2m$ will typically represent a very large sample in
one--dimension.  The estimation in {\tt tgp} uses a moving average
with squared distance weights and a window containing the {\tt
  span}$*2m$ nearest points (the {\tt span} argument defaults to 0.3).


\subsection{Examples}

We illustrate the capabilities of the {\tt sens} function by looking
at the Friedman function considered earlier in this vignette. 
 The function that describes the
responses ($Z$), observed with standard Normal noise, has mean
\begin{equation}
E(Z|\mb{x}) =  10 \sin(\pi x_1 x_2) + 20(x_3 - 0.5)^2 + 10x_4 + 5 x_5.
\label{eq:f1}
\end{equation}
A sensitivity analysis can be based upon any of the available
regression models (e.g., {\tt btlm}, {\tt bgp}, or {\tt btgp}); we
choose to specify {\tt model=btgpllm} for this example.  
The size of each LHS used in the integration scheme is specified
through {\tt nn.lhs}, such that this is equivalent to $m$ in the above
algorithm description.  Thus the number of locations used for
prediction---the size of the random {\tt XX} prediction matrix---is
{\tt nn.lhs*(ncol(X)+2)}.  In addition, the window for moving average
estimation of the main effects is {\tt span*2*nn.lhs} (independent of
this, an {\tt ngrid} argument with a default setting of 
{\tt ngrid=100} dictates the number of grid points in each input
dimension upon which main effects will be estimated).
<<>>=
f <- friedman.1.data(250) 
@ 
This function actually generates 10 covariates, the last five of
which are completely un-influential.  We'll include one of these
($x_6$) to show what the sensitivity analysis looks like for unrelated
variables.
<<>>=
Xf <- f[, 1:6] 
Zf <- f$Y 
sf <- sens(X=Xf, Z=Zf, nn.lhs=600, model=bgpllm, verb=0)
@ 

The progress indicators printed to the screen (for {\tt verb > 0}) are
the same as would be obtained under the specified regression {\tt
  model}---{\tt bgpllm} in this case---so we suppress them here.  All
of the same options (e.g., {\tt BTE}, {\tt R}, etc.) apply, although
if using the {\tt trace} capabilities one should be aware that the
{\tt XX} matrix is changing throughout the MCMC.  The {\tt sens}
function returns a \verb!"tgp"!-class object, and all of the SA related
material is included in the {\tt sens} list within this object.
<<>>=
names(sf$sens)
@ 
The object provides the SA parameters ({\tt par}), the grid of
locations for main effect prediction ({\tt Xgrid}), the mean and
interval estimates for these main effects ({\tt ZZ.mean}, {\tt ZZ.q1},
and {\tt ZZ.q2}), and full posterior via samples of the sensitivity
indices ({\tt S} and {\tt T}).

The plot function for \verb!"tgp"!-class objects now provides a
variety of ways to visualize the results of a sensitivity analysis.
This capability is accessed by specifying {\tt layout="sens"} in the
standard {\tt plot} command.  By default, the mean posterior main
effects are plotted next to boxplot summaries of the posterior sample
for each $S_j$ and $T_j$ index, as in Figure \ref{full}.
\begin{figure}[ht!]
<<label=sens-full,fig=TRUE,echo=TRUE,width=10,height=5,include=FALSE>>= 
plot(sf, layout="sens", legendloc="topleft")
@
<<echo=false,results=hide>>=
graphics.off()
@
\includegraphics[width=6.5in,trim=0 10 0 10]{tgp2-sens-full}
\caption{Full sensitivity analysis results for the Friedman function.}
\label{full}
\end{figure}

 
A further note on the role played by {\tt nn.lhs}: As always, the
quality of the regression model estimate depends on the length of the
MCMC.  But now, the quality of sensitivity analysis is directly
influenced by the size of the LHS used for integral approximation; as
with any Monte Carlo integration scheme, the sample size (i.e., {\tt
  nn.lhs}) must increase with the dimensionality of the problem.  In
particular, it can be seen in the estimation procedure described above
that the total sensitivity indices (the $T_j$'s) are not forced to be
non-negative. If negative values occur it is necessary to increase
{\tt nn.lhs}.  In any case, the {\tt plot.tgp} function 
changes any of the negative values to zero for purposes of illustration.

The {\tt maineff} argument can be used to plot either selected main effects
(Figure \ref{mains}),
\begin{figure}[ht!]
<<label=sens-mains,fig=TRUE,echo=TRUE,width=10,height=4,include=FALSE>>= 
par(mar=c(4,2,4,2), mfrow=c(2,3))
plot(sf, layout="sens", maineff=t(1:6))
@
<<echo=false,results=hide>>=
graphics.off()
@
\centering
\includegraphics[width=6.6in]{tgp2-sens-mains}
\caption{Friedman function main effects, with posterior 90\% intervals.}
  \label{mains}
\end{figure}
or just the sensitivity indices (Figure \ref{indices}).
\begin{figure}[ht!]
<<label=sens-indices,fig=TRUE,echo=TRUE,width=10,height=5,include=FALSE>>= 
plot(sf, layout="sens", maineff=FALSE)
@
<<echo=false,results=hide>>=
graphics.off()
@
\centering
\includegraphics[width=6.5in,trim=0 15 0 15]{tgp2-sens-indices}
\caption{Sensitivity indices for the Friedman function.}
  \label{indices}
\end{figure}
Note that the posterior intervals shown in these plots represent
uncertainty about both the function response and the integration
estimates; this full quantification of uncertainty is not presently
available in any alternative SA procedures.  These plots may be
compared to what we know about the Friedman function (refer to
Eq.~(\ref{eq:f1})) to evaluate the analysis.  The main effects
correspond to what we would expect: sine waves for $x_1$ and $x_2$, a
parabola for $x_3$, and linear effects for $x_4$ and $x_5$.  The
sensitivity indices show $x_1$ and $x_2$ contributing roughly
equivalent amounts of variation, while $x_4$ is relatively more
influential than $x_5$.  Full effect sensitivity indices for $x_3$,
$x_4$, and $x_5$ are roughly the same as the first order indices, but
(due to the interaction in the Friedman function) the sensitivity
indices for the total effect of $x_1$ and $x_2$ are significantly
larger than the corresponding first order indices.  Finally, our SA is
able to determine that $x_6$ is unrelated to the response.

This analysis assumes the default uncertainty
distribution, which is uniform over the range of  input data.  In
other scenarios, it is useful to specify an informative $u(\mb{x})$.
In the {\tt sens} function, properties of $u$ are defined through the
{\tt rect}, {\tt shape}, and {\tt mode} arguments.  To guarantee
integrability of our indices, we have restricted ourselves to bounded
uncertainty distributions.  Hence, {\tt rect} defines these bounds. In
particular, this defines the domain from which the LHSs are to be
taken. We then use independent scaled beta distributions,
parameterized by the {\tt shape} parameter and distribution {\tt
  mode}, to define an informative uncertainty distribution over this
domain.

As an example of sensitivity analysis under an informative uncertainty
distribution, consider the {\tt airquality} data available with the
base distribution of {\sf R}.  This data set contains daily readings
of mean ozone in parts per billion ({\it Ozone}), solar radiation
({\it Solar.R}), wind speed ({\it Wind}), and maximum temperature
({\it Temp}) for New York City, between May 1 and September 30, 1973.
Suppose that we are interested in the sensitivity of air quality to
natural changes in {\it Solar.R},{\it Wind}, and {\it Temp}.  For
convenience, we will build our uncertainty distribution while assuming
independence between these inputs.  Hence, for each
variable, the input uncertainty distribution will be a scaled beta
with {\tt shape=2},
and {\tt mode} equal to the data mean.
<<>>=
X <- airquality[,2:4]
Z <- airquality$Ozone
rect <- t(apply(X, 2, range, na.rm=TRUE))
mode <- apply(X , 2, mean, na.rm=TRUE)
shape <- rep(2,3)
@ 
LHS samples from the uncertainty distribution are shown in Figure
(\ref{udraw})
\begin{figure}[ht!]
<<label=sens-udraw,fig=TRUE,echo=TRUE,width=8,height=4,include=FALSE>>= 
Udraw <- lhs(300, rect=rect, mode=mode, shape=shape)
par(mfrow=c(1,3), mar=c(4,2,4,2))
for(i in 1:3){
  hist(Udraw[,i], breaks=10,xlab=names(X)[i], 
       main="",ylab="", border=grey(.9), col=8) 
}  
@
<<echo=false,results=hide>>=
graphics.off()
@
\centering
\includegraphics[width=6in,trim=0 0 0 30]{tgp2-sens-udraw}
\caption{A sample from the marginal uncertainty distribution for the
  airquality data.}
\label{udraw}
\end{figure}

Due to missing data (discarded in the current version of {\tt tgp}),
we suppress warnings for the sensitivity analysis.  We shall use
the default {\tt model=btgp}.
<<>>=
s.air <- suppressWarnings(sens(X=X, Z=Z, nn.lhs=300, rect=rect, 
                               shape=shape, mode=mode, verb=0))
@ 
Figure (\ref{air1}) shows the results from this analysis.
\begin{figure}[ht!]
<<label=sens-air1,fig=TRUE,echo=TRUE,width=10,height=5,include=FALSE>>= 
plot(s.air, layout="sens")
@
<<echo=false,results=hide>>=
graphics.off()
@
\centering
\includegraphics[width=6.5in,trim=0 15 0 15]{tgp2-sens-air1}
\caption{Sensitivity of NYC airquality to natural variation in wind,
  sun, and temperature.}
\label{air1}
\end{figure}


Through use of {\tt predict.tgp}, it is possible to quickly
re-analyze with respect to a new uncertainty distribution without
running new MCMC.  We can, for example, look at sensitivity for air
quality on only low--wind days.  We thus alter the uncertainty
distribution (assuming that things remain the same for the other
variables)
<<>>=
rect[2,] <- c(0,5)
mode[2] <- 2
shape[2] <- 2
@ 
and build a set of parameters {\tt sens.p} with the {\tt sens}
function by setting {\tt model=NULL}.  
<<>>=
sens.p <- suppressWarnings(sens(X=X,Z=Z,nn.lhs=300, model=NULL, rect=rect, shape=shape, mode=mode))
@ 
\begin{figure}[ht!]
<<label=sens-air2,fig=TRUE,echo=TRUE,width=10,height=5,include=FALSE>>= 
s.air2 <- predict(s.air, BTE=c(1,1000,1), sens.p=sens.p, verb=0) 
plot(s.air2, layout="sens")
@
<<echo=false,results=hide>>=
graphics.off()
@
\includegraphics[width=6.5in,trim=0 15 0 15]{tgp2-sens-air2}
\caption{Air quality sensitivity on low-wind days.}
\label{air2}
\end{figure}
Figures (\ref{air1}) and (\ref{air2}) both show total effect indices
which are much larger than the respective first order sensitivities.
As one would expect, the effect on airquality is manifest largely
through an interaction between variables.

Finally, it is also possible to perform SA with binary covariates,
included in the regression model as described in Section 1.  In this
case, the uncertainty distribution is naturally characterized by a
Bernoulli density. Setting {\tt shape[i]=0} informs {\tt sens} that
the relevant variable is binary (perhaps encoding a categorical input
as in Section \ref{sec:cat}), and that the Bernoulli uncertainty
distribution should be used.  In this case, the {\tt mode[i]}
parameter dictates the probability parameter for the Bernoulli, and we
must have {\tt rect[i,] = c(0,1)}.  As an example, we re-analyze the
original air quality data with temperature included as an indicator
variable (set to one if temperature > 79, the median, and zero
otherwise).
<<>>=
X$Temp[X$Temp >70] <- 1
X$Temp[X$Temp >1] <- 0
rect <- t(apply(X, 2, range, na.rm=TRUE))
mode <- apply(X , 2, mean, na.rm=TRUE)
shape <- c(2,2,0)
s.air <- suppressWarnings(sens(X=X, Z=Z, nn.lhs=300, rect=rect, 
                               shape=shape, mode=mode, verb=0, basemax=2))
@ 
\begin{figure}[ht!]
<<label=sens-air3,fig=TRUE,echo=TRUE,width=10,height=5,include=FALSE>>= 
plot(s.air, layout="sens")
@
<<echo=false,results=hide>>=
graphics.off()
@
\centering
\includegraphics[width=6.5in,trim=0 15 0 15]{tgp2-sens-air3}
\caption{Sensitivity of NYC airquality to natural variation in wind,
  sun, and a binary temperature variable (for a threshold of 79 degrees).}
\label{air3}
\end{figure}
Figure (\ref{air3}) shows the results from this analysis.


