scar: Compute the maximum likelihood estimator of the generalised additive regression with shape constraints

Description

This function uses the active set algorithm to compute the maximum likelihood estimator (mle) of the generalised additive regression with shape constraints. Each component function of the additive predictors is assumed to belong to one of the nine possible shape restrictions. The estimator's value at the data points is unique.

The output is an object of class scar which contains all the information needed to plot the estimator using the plot method, or to evaluate it using the predict method.

Usage

scar(x, y, shape = rep("l", d), family = gaussian(),
  weights = rep(1, length(y)), epsilon = 1e-08)

Arguments

Observed covariates in $R^{d}$ , in the form of an $n \times d$ numeric matrix.

Observed responses, in the form of a numeric vector of length $n$ .

shape

A vector that specifies the shape restrictions for each component function, in the form of a string vector of length $d$ . The string allowed and its corresponding shape constraint is listed as follows (see Details):

l: linear

in: monotonically increasing

de: monotonically decreasing

cvx: convex

cvxin: convex and increasing

cvxde: convex and decreasing

ccv: concave

ccvin: concave and increasing

ccvde: concave and decreasing

family

A description of the error distribution and link function to be used in the model. This can be a character string naming a family function, a family function or the result of a call to a family function. Currently only the following five common exponential families are allowed: Gaussian, Binomial, Poisson, and Gamma. By default the canonical link function is used.

weights

An optional vector of prior weights to be used when maximising the likelihood. It is a numeric vector of length $n$ . By default equal weights are used.

epsilon

Positive convergence tolerance epsilon when performing the iteratively reweighted least squares (IRLS) method at each iteration of the active set algorithm. See glm.control for more details.

Value

An object of class scar, with the following components:

Covariates copied from input.

Response copied from input.

shape

Shape vector copied from input.

weights

The vector of weights copied from input.

family

The exponential family copied from input.

componentfit

Value of the fitted component function at each observed covariate, in the form of an $n \times d$ numeric matrix, where the element at the $i$ -th row and the $j$ -th column is the value of $f_{j}$ at the $j$ -th coordinate of $X_{i}$ , with the identifiability condition satisfied (see Details)

constant

The estimated value of the constant $c$ in the additive function $f$ (see Details).

deviance

Up to a constant, minus twice the maximised log-likelihood. Where applicable, the constant is chosen to make the saturated model to have zero deviance. See also glm.

nulldeviance

The deviance for the null model.

iter

Total number of iterations of the active set algorithm

Details

For $i = 1, \dots, n$ , let $X_{i}$ be the $d$ -dimensional covariates, $Y_{i}$ be the corresponding one-dimensional response and $w_{i}$ be its weight. The generalised additive model can be written as $g (μ) = f (x),$ where $x = (x_{1}, \dots, x_{d})^{T}$ , $g$ is a known link function and $f$ is an additive function (to be estimated).

Assume the canonical link function is used here, then the maximum likelihood estimator of the generalised additive model based on observations $(X_{1}, Y_{1}), \dots, (X_{n}, Y_{n})$ is the function that maximises $\frac{1}{n} \sum_{i = 1}^{n} w_{i} {Y_{i} f (X_{i}) - B (f (X_{i}))}$ subject to the restrictions that for every $j = 1, \dots, d$ , the $j$ -th additive component of $f$ satisfies the constraint indicated by the $j$ -th element of shape. Here $B (.)$ is the log-partition function of the specified exponential family distribution, and $w_{i}$ are the weights. For i.i.d. data, $w_{i}$ should be $1$ for each $i$ .

To make each component of $f$ identifiable, we write $f (x) = \sum_{j = 1}^{d} f_{j} (x_{j}) + c$ and let $f_{j} (0) = 0$ for every $j = 1, \dots, d$ . In case zero is outside the range of the $j$ -th observed covariate, for the sake of convenience, we set $f_{j}$ to be zero at the sample mean of the $j$ -th predictor.

This problem can then be re-written as a concave optimisation problem, and our function uses the active set algorithm to find out the maximum likelihood estimator. A general introduction can be found in Nocedal and Wright (2006). A detailed description of our algorithm can be found in Chen and Samworth (2016). See also Groeneboom, Jongbloed and Wellner (2008) for some theoretical supports.

References

Chen, Y. and Samworth, R. J. (2016). Generalized additive and index models with shape constraints. Journal of the Royal Statistical Society: Series B, 78, 729-754.

Groeneboom, P., Jongbloed, G. and Wellner, J.A. (2008). The support reduction algorithm for computing non-parametric function estimates in mixture models. Scandinavian Journal of Statistics, 35, 385-399.

Hastie, T. and Tibshirani, R. (1990). Generalized Additive Models. Chapman and Hall, London.

Meyer, M. C. (2013). Semi-parametric additive constrained regression. Journal of nonparametric statistics, 25, 715-743.

Nocedal, J., and Wright, S. J. (2006). Numerical Optimization, 2nd edition. Springer, New York.

Robertson, T., Wright, F. T. and Dykstra, R. L. (1988). Order Restricted Statistical Inference. Wiley, New York.

Venables, W. N. and Ripley, B. D. (2002). Modern Applied Statistics with S. Springer, New York.

Wood, S.N. (2004). Stable and efficient multiple smoothing parameter estimation for generalized additive models. Journal of American Statistical Association, 99, 673-686.

Examples

Run this code

# NOT RUN {
## An example in the Poission additive regression setting:
## Define the additive function f on the scale of the predictors
f<-function(x){
  return(1*abs(x[,1]) + 2*(x[,2])^2 + 3*(abs(x[,3]))^3) 
}

## Simulate the covariates and the responses
## covariates are drawn uniformly from [-1,1]^3
set.seed(0)
d = 3
n = 500
x = matrix(runif(n*d)*2-1,nrow=n,ncol=d) 
rpoisson <- function(m){rpois(1,exp(m))}
y = sapply(f(x),rpoisson)

## All the components are convex so one can use scar
shape=c("cvx","cvx","cvx")
object = scar(x,y,shape=shape, family=poisson())

## Plot each component of the estimatied additive function
plot(object)

## Evaluate the estimated additive function at 10^4 random points 
## drawing from the interior of the support
testx = matrix((runif(10000*d)*1.96-0.98),ncol=d)
testf = predict(object,testx)

## and calculate the (estimated) absolute prediction error
mean(abs(testf-f(testx))) 
# }

Run the code above in your browser using DataLab

Get 50% off unlimited learning