SSM: A package for fitting smooth supersaturated models (SSM).

Description

The SSM package provides an S4 class to represent smooth supersaturated models, along with functions for fitting SSM to data, computing Sobol indices for the SSM and estimating metamodel error with a Gaussian process.

SSM functions

There are three important functions in the package.

fit.ssm returns an SSM object that fits an SSM to the data and, by default, computes the Sobol indices, and Total interaction indices of the model. Optionally, the metamodel error can be estimated using a Gaussian process. The fitted SSM smooths over [-1, 1]^d and uses a hierarchical basis of 20*d+n Legendre polynomials by default but the function is highly customisable.

predict.SSM returns the model prediction at a point, including a credible interval if a metamodel error GP has been fit.

plot.SSM plots the main effects of the SSM.

An S4 class to represent a smooth supersaturated model

Description

An S4 class to represent a smooth supersaturated model

Slots

dimension

A number indicating the number of variables in the design.

design

A matrix with rows indicating the design point and columns indicating the variable.

design_size

A number indicating the number of design points.

response

A design_size length vector of responses.

theta

A vector containing the fitted model coefficients.

basis

A matrix with each row being the exponent vector of a polynomial term.

basis_size

A number indicating the number of basis terms used in the model. This may be different from nrow(basis) if terms are excluded.

include

A vector containing the row numbers of the basis polynomials used in the model. This is used when interactions or variables are being excluded from the model.

K

A semi-positive definite matrix that defines the smoothing criteria.

P

A matrix that defines the polynomial basis in terms of a monomial basis.

design_model_matrix

A matrix.

variances

A vector of length basis_size containing the term variances.

total_variance

A length one vector containing the total variance.

main_sobol

A dimension length vector containing the Sobol index for each variable.

main_ind

A logical matrix indicating whether each term is included in the main effect corresponding to the column.

total_sobol

A dimension length vector containing the Total sensitivity index for each variable.

total_ind

A logical matrix indicating whether each term is included in the Total sensitivity index corresponding to the column.

int_sobol

A vector containing the Sobol index for interactions.

int_factors

A list of length the same as int_sobol indicating which interaction corresponds with each entry in int_sobol.

total_int

A vector containing the Total interaction indices of all second order interactions.

total_int_factors

A matrix where each row indicates the variables associated with the corresponding interaction in total_int.

distance

A matrix containing the distances used for computing the covariance matrix of the GP metamodel error estimate.

distance_type

A character defining the distance type used for computing distance. Can be one of "distance", "line", "product", "area", "proddiff", or "smoothdiff".

type

A character, either "exp", "matern32", that selects the correlation function used for the GP metamodel error estimate.

covariance

A positive definite matrix. The covariance matrix of the GP metamodel error estimate prior to scaling by sigma.

residuals

A design_size length vector containing the Leave-One-Out errors of the model at each design point.

sigma

A number indicating the scaling factor for covariance.

r

A number indicating the length factor for the correlation function.

local_smoothness

A design_size length vector containing the model smoothness at each design point.

LOO_RMSE

A number. The Leave-One-Out root mean square error.

legendre

logical. Indicates whether the default Legendre polynomial basis is being used.

fail

logical. Indicates whether the model fit was successful.

Generate all desired exponent vectors of a given degree.

Description

comb recursively generates all the desired exponent vectors of a given degree and is called by degl to generate the matrix put into the basis slot of the SSM object.

Usage

comb(d, deg, N = choose(d + deg - 1, deg), vec, start = TRUE, parent,
  exclude = list())

Arguments

Details

This function is called by degl during the process of constructing the objects related to the basis of a smooth supersaturated model. It constructs the first N exponent vectors of degree d, excluding those which are non-zero only in the columns specified by vectors listed in exclude. It operates recursively.

Value

A matrix of exponent vectors

Compute the unscaled covariance matrix.

Description

Computes the unscaled covariance matrix for an SSM object for a given length parameter.

Usage

compute.covariance(ssm, r)

Arguments

Details

The covariance matrix uses the correlation function specified by the SSM slot type. "matern32" uses a Matern 3/2 correlation while "exp" uses a squared exponential.

Value

The unscaled covariance matrix.

Compute unscaled covariance matrix from a supplied distance matrix and length parameter.

Description

This is a legacy function and is not used.

Usage

compute.covariance.from.distance(distance, r, type = "exp", ssm)

Arguments

Compute Total interaction indices and Sobol indices for higher order interactions.

Description

This function computes the Total interaction indices of all second order interactions present in the model. If d < 11 then it also computes the Sobol indices of all interactions (second order and higher) in the model. It identifies the relevant rows of the SSM basis slot for each interaction and sums the term variances. Called by update.sensitivity.

Usage

compute.interactions(ssm)

Arguments

Value

An SSM object. This is the ssm input with the relevant slots updated.

Compute main effects

Description

This function computes the Sobol indices of the model variables, also known as the main effects. It identifies the relevant rows of the SSM basis slot and sums the term variances. Called by update.sensitivity.

Usage

compute.main.effects(ssm)

Arguments

Value

An SSM object. This is the ssm input with the relevant slots updated.

Compute the Leave-One-Out error at all design points.

Description

This function repeatedly fits the SSM to all but one held out point in turn and computes the prediction error at the held out point. These errors are collated into a vector indicating the Leave-One-Out error at each design point.

Usage

compute.residuals(ssm)

Arguments

Value

A vector of Leave-One_Out errors.

Compute the Sobol index for a given interaction.

Description

This computes the Sobol index for a given interaction. The relevant term variances are identified and summed and the resulting variance normalized and returned.

Usage

compute.specific.interaction(factors, ssm)

Arguments

Value

A number. The Sobol index of the requested interaction.

Compute Total interaction variance

Description

This computes the Total interaction index for a given interaction. The relevant term variances are identified and summed and the resulting variance normalized and returned.

Usage

compute.specific.total.interaction(factors, ssm)

Arguments

Value

A number. The Total interaction index of the requested interaction.

Compute Total effects

Description

This function computes the Total indices of the model variables. It identifies the relevant rows of the SSM basis slot and sums the term variances. Called by update.sensitivity.

Usage

compute.total.effects(ssm)

Arguments

Value

An SSM object. This is the ssm input with the relevant slots updated.

Compute the concentrated likelihood of a covariance matrix.

Description

Computes the concentrated likelihood of the covariance matrix of an SSM object, given a length parameter and the SSM Leave-One-Out errors.

Usage

concentrated.likelihood(r, ssm)

Arguments

Value

A number. The concentrated likelihood.

Construct the K matrix for a given multivariate basis.

Description

This function constructs the K matrix for a given multivariate basis assuming the basis is a Legendre polynomial basis and the smoothing criterion is the Frobenius norm of the Hessian integrated over [-1, 1]^d.

Usage

construct.K(basis)

Arguments

Value

A matrix where each entry is <f, g> with

<f, g> = \int_X \sum_{i,j}\frac{d^2f}{dx_idx_j}\frac{d^2g}{dx_idx_j} dx,

with f, g being the Legendre polynomials described by the appropriate exponent vectors.

Construct the K matrix for a given univariate basis.

Description

This function constructs the K matrix for a given multivariate basis assuming the basis is a Legendre polynomial basis and the smoothing criterion is the Frobenius norm of the Hessian integrated over [-1, 1].

Usage

construct.K.1d(basis)

Arguments

Value

A matrix where each entry is <f, g> with

<f, g> = \int_X \frac{d^2f}{dx^2}\frac{d^2g}{dx^2} dx,

with f, g being the Legendre polynomials described by the appropriate exponent vectors.

Construct the change of basis matrix from multivariate monomials to Legendre polynomials.

Description

Let f(x) be the vector of monomials implied by the matrix basis. Let P be the matrix generated by Legendre_P. Then g(x)=P^Tf(x) is the vector of Legendre polynomials with leading terms corresponding to basis.

Usage

construct.P(basis)

Arguments

Value

A matrix which functions as a change of basis from monomials to Legendre.

Construct the change of basis matrix from univariate monomials to Legendre polynomials.

Description

Let f(x) be the vector of monomials implied by the matrix basis. Let P be the matrix generated by Legendre_P_1d. Then g(x)=P^Tf(x) is the vector of Legendre polynomials with leading terms corresponding to basis. e.g. the matrix (0 1 2)^T implies f(x)^T = (1 x x^2). Then g(x)^T = (L_0, L_1, L_2).

Usage

construct.P.1d(basis)

Arguments

Value

A matrix which functions as a change of basis from monomials to Legendre.

Construct the design model matrix

Description

Constructs the design model matrix corresponding to the given design and basis. Used by fit.ssm.

Usage

construct.dmm(basis, design, P)

Arguments

Details

The argument basis defines a monomial basis and the change of basis matrix P is used to convert this to a polynomial basis. The function returns the design model matrix corresponding to this polynomial basis and the design given by design. Note that if P is the appropriately sized identity matrix then the model basis a monomial one.

Value

The n x N design model matrix.

Construct matrix of exponent vectors.

Description

This function is called by fit.ssm when it needs to construct the matrix of exponent vectors stored in the basis slot.

Usage

degl(d, N, exclude = list())

Arguments

Details

The algorithm works by repeated calls to comb to generate all possible exponent vectors of a given degree until N vectors have been generated. Any generated vector is checked to make sure that it's non-zero entries do not match a vector provided in exclude before being added to the output matrix.

Estimate the parameters of the metamodel error estimating GP.

Description

This function estimates the parameters of the metamodel error estimating Gaussian process using maximum likelihood methods to identify the length parameter r of the given correlation function.

Usage

estimate.GP(ssm, type)

Arguments

Details

Since the concentrated likelihood function is often flat, this function calls optimize several times using different search intervals to avoid cases when the optimize algorithm misses maxima when the interval is too large.

The scaling parameter sigma is found analytically once r has been determined.

Value

An SSM object that is the same as the input except with estimates for r and sigma in the appropriate slots.

Compute the SSM vector of parameters.

Description

This function is called by fit.ssm to compute the coefficient vector that interpolates the data and minimises the smoothness of the resulting model.

Usage

find.theta(response, K, design_model, tol = .Machine$double.eps)

Arguments

Value

A vector of parameters of length N if the model is fit successfully. NA is returned should solve not invert the required matrix.

Fit a smooth supersaturated model

Description

fit.ssm fits a smooth supersaturated model to given data. By default the model smooths over [-1, 1]^d and uses a basis of Legendre polynomials. Many model parameters such as basis size and smoothing criterion can be user-defined. Optionally, sensitivity indices and a metamodel error estimating Gaussian process can be computed.

Usage

fit.ssm(design, response, ssm, basis, basis_size, K, P, design_model_matrix,
  SA = FALSE, GP = FALSE, type = "exp", validation = FALSE,
  exclude = list(), distance_type = "distance")

Arguments

Details

Returns an SSM object containing the fitted smooth supersaturated model. Minimal arguments required are design and response. This will result in a model using Legendre polynomials and smoothing over [-1, 1]^d. All other arguments will be assigned automatically if not specified.

If the model is unable to be fit due to numerical instability a warning will be given and the returned SSM object will have a zero vector of appropriate length in the theta slot. The basis_size parameter is often useful in this situation. Try reducing the basis_size until a model fit is successful. Ideally the basis size should be as large as possible without causing instability in the predictions (see the example below).

If SA is TRUE then sensitivty analysis will be performed on the model. Sobol indices for main effects, Total indices, and Total interaction indices for second order interactions will be computed and held in the slots main_sobol, total_sobol and total_int respectively. If the number of factors is < 11 then Sobol indices will be computed for all order interactions and stored in int_sobol. Default behaviour is to assume each input is uniformly distributed over [-1, 1]. If P has been used-defined then the polynomials defined by P are assumed to be an orthonormal system with respect to some measure on the reals. See update.sensitivity for more details.

If GP is TRUE (default behaviour is false due to computational cost) then a the metamodel error is estimated using a zero-mean Gaussian process with a constant trend. Scaling and length parameters are estimated using maximum likelihood methods using the Leave-One-Out model errors and stored in the sigma and r slots respectively. Model predictions using predict.SSM will then include credible intervals. The distance between points is defined by the distance_type argument and is set to the standard Euclidean distance by default. See new.distance for other options although they are experimental and subject to erratic behaviour. The default correlation function used is the square exponential although this can be changed to a Matern 3/2 function by setting the type argument to "matern32".

If validation is TRUE then the Leave-One_Out error at each design point will be computed and stored in the residuals slot, and the LOO RMSE computed and stored in the LOO_RMSE slot. Note that if GP is TRUE then these values will be computed regardless of the value of validation as they are required to fit the metamodel error estimate GP.

Value

An SSM object.

See Also

predict.SSM for model predictions for SSM, and plot.SSM for plotting main effects of SSM. transform11 is useful for transforming data to [-1, 1]^d.

Examples

# A simple one factor example
X <- seq(-1,1,0.5) # design
Y <- c(0,1,0,0.5,0) # response
s <- fit.ssm(X,Y)
s
plot(s)
predict(s,0.3)

# used defined basis sizes

# A model that is too large to fit
## Not run: 
s <- fit.ssm(X, Y, basis_size=80)

## End(Not run)
# A large model that can be fit but is unstable
s <- fit.ssm(X, Y, basis_size=70)
plot(s)
# A model larger than default that is not unstable
s <- fit.ssm(X, Y, basis_size=40)
plot(s)

# with metamodel error estimate

s <- fit.ssm(X, Y, GP=TRUE)
plot(s)
predict(s,0.3)

# Sensitivity analysis and main effect plotting

# A design of 20 points over [-1, 1]^d
X <- matrix(runif(20, -1, 1), ncol = 2)
Y <- runif(10)
s <- fit.ssm(X, Y, SA = TRUE)
s
sensitivity.plot(s)
plot(s)

Compute entry of K matrix.

Description

Called by construct.K to compute entries of the K matrix.

Usage

get.K.element(n, m, basis)

Arguments

Value

The inner product of the n, m-th basis elements as defined for construct.K.

Identify main effect terms

Description

This is used by the sensitivity computation functions to identify which basis terms are associated with each main effect.

Usage

## S3 method for class 'main.effect.terms'
identify(i, basis)

Arguments

Value

A logical vector. Each term corresponds with the associated row of basis.

Identify total effect terms

Description

This is used by the sensitivity computation functions to identify which basis terms are associated with each total effect.

Usage

## S3 method for class 'total.effect.terms'
identify(i, basis)

Arguments

Value

A logical vector. Each term corresponds with the associated row of basis.

Plot the concentrated likelihood of an SSM.

Description

Plot the concentrated likelihood used to estimate the parameters of the metamodel error estimating Gaussian process.

Usage

likelihood.plot(ssm, xrange = c(0, 1000), grid = 200)

Arguments

Details

As a diagnostic it can be helpful to look at the concentrated likelihood function of the correlation function used to estimate the metamodel error. Flat likelihood functions make it difficult to pick a suitable r length parameter. Note that r and sigma can be set manually.

Examples

data(attitude)
X <- transform11(attitude[ 2:7])
Y <- attitude[ , 1]
s <- fit.ssm(X, Y, GP = TRUE)
likelihood.plot(s)
likelihood.plot(s, xrange = c(0, 20))

Average the values in a vector between two cutoff points specified by a separate vector.

Description

The values of smoothnessdesign that are in the interval specifed by x are identified and the associated values in smoothness are averaged. This is used in the computation of line integrals of smoothness used by a particular distance measure when computing metamodel error estimates.

Usage

lineij(x, smoothness, smoothnessdesign)

Arguments

Value

A number.

Compute the distance matrix of an SSM design.

Description

Computes the distance matrix associated with the design held in the design slot of an SSM object. Used in the construction of the metamodel error estimating Gaussian process.

Usage

new.distance(type = "distance", ssm, line.grid = 100)

Arguments

Details

Implemented types of distance measure are:

• "distance" Standard Euclidean distance.

• "line" The line integral of the smoothness between design points. This is not implemented for data where d > 1.

• "product" The product of the Euclidean distance and the local smoothness at both points.

• "area" The sum of the local smoothness at the points multiplied by the Euclidean distance.

• "proddiff" Multiplies the difference in local smoothness between points by the Euclidean distance.

• "smoothdiff" The difference between the local smoothness of points.

All measures other than "distance" are experimental and should be used with caution.

Value

A matrix.

Optimize concentrated likelihood.

Description

Wrapper for optimizing concentrated.likelihood. Used by estimate.GP.

Usage

optimize.by.interval.maximum(int, ...)

Arguments

Compute second partial derivative of a smooth supersaturated model at all design points.

Description

Computes a second partial derivative of a smooth supersaturated model at a design point. Used in the computation of distance measures based on local smoothness.

Usage

partial.deriv.ssm(indices, ssm)

Arguments

Value

A single column vector containing the second partial derivatives with respect to the requested variables, evaluated at all design points.

Plot smooth supersaturated model main effects

Description

plot.SSM is a plot method for SSM objects. It plots the main effects of the SSM only, that is the subset of basis terms that are dependent on a single variable only. For single variable data this is a plot of the complete model.

Usage

## S3 method for class 'SSM'
plot(x, ..., grid = 200, yrange = "full", GP = TRUE)

Arguments

Details

For each variable, the effect is plotted over [-1, 1] by default although passing an alternate range to the xlim argument will override this.

The yrange argument is designed to automatically compute the relevant plot range for each effect. By default a ylim value is passed to plot that covers the range of responses. "full" results in a ylim value that covers the range of predictions or, if appropriate, the range of the metamodel error credible interval.

For single variable data, setting GP to TRUE will plot a credible interval for the metamodel error estimating Gaussian process if this has been computed for the SSM object.

Examples

# A single variable example
X <- seq(-1, 1, 0.25)
Y <- sapply(X, "^", 3)
s <- fit.ssm(X, Y, GP = TRUE)
plot(s)

# A six variable example
data(attitude)
X <- transform11(attitude[ 2:7])
Y <- attitude[ , 1]
s <- fit.ssm(X, Y)
plot(s)

Point prediction of smooth supersaturated models.

Description

This method gives the prediction of an SSM object at a point. If the SSM has a metamodel error estimate then a (1 - \alpha) credible interval is also output.

Usage

## S3 method for class 'SSM'
predict(object, x, alpha = 0.05, ...)

Arguments

Value

Either a number if the SSM has no metamodel error estimating Gaussian process, or three numbers giving the model prediction ($model), and the lower and upper bounds of the credible interval ($lower and $upper) respectively.

Examples

data(attitude)
X <- transform11(attitude[ 2:7])
Y <- attitude[ , 1]
# with no metamodel error estimating GP.
s <- fit.ssm(X, Y)
predict(s, rep(1,6))

# with metamodel error estimating GP.
s <- fit.ssm(X, Y, GP = TRUE)
predict(s, rep(1,6))

Plot the sensitivity indices of a smooth supersaturated model.

Description

sensitivity.plot visualises the sensitivity indices of a given smooth supersaturated model using barplot. If the SA flag was not set to TRUE when fit.ssm was run to fit the model, then update.sensitivity should be used to compute the model variances. If not, this function will return an error message.

Usage

sensitivity.plot(ssm, type = "main_total", ...)

Arguments

Details

There are four classes of plot available:

• "sobol" Produces a barplot of all Sobol indices. If there are more than 10 factors then Sobol indices will not have been computed for interactions and only the Sobol indices for main effects will be plotted. Main effects are in red, interactions are in grey.

• "main_sobol" Produces a barplot of Sobol indices for main effects only.

• "main_total" Produces a barplot of Total indices for main effects only. This is the default behaviour.

• "total" Produces a barplot of Total indices for main effects and all second order interactions. Main effects are in red, interactions are in grey.

Note that variables and interactions are not labelled in the plots since there can be too many bars to label clearly.

Examples

# A 20 point design in four variables
X <- matrix(runif(80, -1, 1), ncol = 4)
Y <- runif(20)
s <- fit.ssm(X, Y, SA = TRUE)
sensitivity.plot(s)

# In the next plots, the grey bars indicate interactions.
sensitivity.plot(s, "sobol")
sensitivity.plot(s, "total")
# Identifying particular indices is best done using the information held in
# the SSM object.  The following orders s@total_int_factors so the
# interaction indicated by the top row is the most important.
ind <- order(s@total_int, decreasing = TRUE)
s@total_int_factors[ind, ]

Summarise SSM class object

Description

Printing an SSM will summarise the parameters d, N, and n. If no model has been fit then this will be noted, otherwise The smoothness of the SSM will be shown. If sensitivity analysis has been performed, the Sobol indices and Total indices for main effects are displayed and if cross-validation has been performed the Standardised LOO RMSE is also shown.

Usage

## S4 method for signature 'SSM'
show(object)

Arguments

Compute the smoothness of an SSM at all design points.

Description

This function evaluates the Frobenius norm of the Hessian matrix at all design points. Used in the compuation of distance measures based on local smoothness.

Usage

smoothness.over.design(ssm)

Arguments

Value

A vector of integers containing the smoothness at each design point.

Transform a design to [-1, 1]^d

Description

This function transforms a design (supplied as a matrix) into the space [-1, 1]^d. This has numerical and computational advantages when using smooth supersaturated models and is assumed by the default fit.ssm behaviour.

Usage

transform11(design)

Arguments

Value

A matrix where each column contains values in [-1, 1]^d.

Examples

X <- transform11(quakes[, 1:4])
apply(X, 2, range)

Update an SSM object with the term variances and Sobol indices

Description

This function computes the term variances, Sobol indices, Total indices and Total interaction indices of a given SSM class object.

Usage

## S3 method for class 'sensitivity'
update(ssm)

Arguments

Details

This function has two modes. If the legendre slot in the SSM object is TRUE (i.e. the default P matrix has been used to fit the SSM) then all the sensitivity indices are computed assuming inputs are uniformly distributed over [-1, 1]^d. If legendre is FALSE (i.e. a user-defined P has been supplied) then the polynomials defined by P are orthonormal with respect to some probability measure and the sensitivity indices are computed assuming that measure defines the distribution of the inputs.

The returned SSM object has term variances held in the variances slot, ordered to match the exponent vectors in the basis slot. The Sobol indices are placed in main_sobol and the Total indices in total_sobol. The Total interaction indices are placed in the total_int slot with identifying labels stored in total_int_factors. If there are less than eleven variables then interaction indices are computed for all order interactions and stored in int_sobol with identifying labels in int_factors.

Value

An SSM object.