Shape-Constrained Additive (index) Regression: a maximum likelihood approach

Description

scar computes the maximum likelihood estimator of the generalised additive and index regression with shape constraints. Each component of the additive function of the predictors is assumed to belong to one of the nine possible shape restrictions: linear, increasing, decreasing, convex, convex and increasing, convex and decreasing, concave, concave and increasing, or concave and decreasing. For the generalised additive regression, the problem is transformed into a convex optimisation problem and the active set algorithm is used to find the optimum. We emphasise that unlike most of the other nonparametric methods, this approach is free of tuning parameters.

Furthermore, we can extend our findings to the generalised additive index regression, where a stochastic search algorithm is proposed to solve the problem.

Details

This package contains a selection of functions for maximum likelihood estimation of the generalised additive (and additive index) regression under shape constraints:

scar computes the maximum likelihood estimator (specified via its value at the observed covariates). Output is a list of class scar which is used as input to various auxiliary functions.

plot.scar produces plots of the maximum likelihood estimator produced by scar on the scale of the additive predictors.

predict.scar obtains predictions either on the scale of the additive predictors or on the scale of the response variable from a fitted scar object.

scair tries to find the maximum likelihood estimator (specified via its value at the observed indices). Output is a list of class scair which is used as input to various auxiliary functions.

plot.scair produces plots of the maximum likelihood estimator produced by scair on the scale of the additive index predictors.

predict.scair obtains predictions either on the scale of the additive index predictors or on the scale of the response variable from a fitted scair object.

The methods proposed here were applied to the following datasets: PhDPublications, decathlon.

Note

The authors gratefully acknowledge the assistance of Ming Yuan for his insights into this problem.

Thanks also go to Mary Meyer for kindly providing the authors with her manuscript prior to its publication.

Author(s)

Yining Chen (maintainer) y.chen101@lse.ac.uk

Richard Samworth r.samworth@statslab.cam.ac.uk

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.

See Also

scar, scair, scam, glm

Examples

## See examples provided in functions scar and scair

Doctoral Publications

Description

Data on the scientific productivity of PhD students in biochemistry.

Usage

data("PhDPublications")

Format

A data frame containing 915 observations on 6 variables:

articles

Number of articles published during last 3 years of PhD.

gender

Factor indicating gender.

married

Factor indicating the marital status.

kids

Number of children less than 6 years old.

prestige

Prestige of the graduate program.

mentor

Number of articles published by student's mentor.

Source

Imported directly from CRAN package AER (Kleiber and Zeileis, 2013) for the sake of convenience.

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.

Kleiber, C. and Zeileis, A. (2013). AER: applied econometrics with R. R package version 1.2-0, http://cran.r-project.org/web/packages/AER/

Long, J.S. (1990). Regression Models for Categorical and Limited Dependent Variables. Thousand Oaks: Sage Publications.

Long, J.S. (1997). The origin of sex differences in science. Social Forces, 68, 1297–1315.

Examples

## Some data pre-processing
data(PhDPublications)
dat = matrix(0,ncol=4,nrow=nrow(PhDPublications))
dat[,1] = as.numeric(PhDPublications$gender == "female")
dat[,2] = as.numeric(PhDPublications$married == "yes")
dat[,3] = as.numeric(PhDPublications$kids)
dat[,4] = as.numeric(PhDPublications$mentor)
y = PhDPublications$articles

## Run scar on the dataset
object = scar(dat,y,shape=c("l","l","de","ccv"),family=poisson())

## Check the effects of mentor
plot(object,which=c(4))

Men's decathlon athletes in 2012

Description

This dataset consists of men's decathlon athletes who scored at least 6500 points in at least one athletic competition in 2012 and scored points in every event there. To avoid data dependency, we include only one performance from each athlete, namely their 2012 personal best (over the whole decathlon).

Usage

data("decathlon")

Format

A data frame containing 614 observations on 10 variables. Different rows represent results of different athletes.

X100m

Points scored in 100 metres.

LJ

Points scored in long jump.

SP

Points scored in shot put.

HJ

Points scored in high jump.

X400m

Points scored in 400 metres.

X110H

Points scored in 110 metres hurdles.

DT

Points scored in discus throw.

PV

Points scored in polt vault.

JT

Points scored in javelin throw.

X1500m

Points scored in 1500 metres.

The point system (on how to convert results to scores) can be found at

http://en.wikipedia.org/wiki/Decathlon

Source

Compiled from the dataset decathlon_raw.

Examples

data(decathlon)
dat = as.matrix(decathlon[,c("X100m","SP","DT")]); y = decathlon$JT

## We consider the problem of predicting a decathlete's javelin performance 
## from their performances in the other decathlon disciplines.
## As an illustration, we only select 100m, shot put and discus throw
## also, we reduce the number of iterations here to shorten the run time
set.seed(1)
object = scair(dat,y,shape=c("in","in"),iter=40, allnonneg=TRUE)
object$index

## Here we can treat the obtained index matrix as a warm start and perform 
## further optimization using e.g. the default R optimisation routine.
## The result is actually only slightly different from the warm start.
## The optimisation itself takes around 5-10 seconds.
## Code provided but commented out here because CRAN team requires all the 
## examples running here to be completed in a few seconds.
# fn<-function(w){
#  wnew = matrix(0,nrow=3,ncol=2)
#  wnew[1,1] = w[1]
#  wnew[2,1] = w[2]
#  wnew[3,1] = 1 - abs(w[1]) - abs(w[2]) 
#  wnew[1,2] = w[3]
#  wnew[2,2] = w[4]
#  wnew[3,2] = 1 - abs(w[3]) - abs(w[4])
#  if ((wnew[3,1] < 0)|| (wnew[3,2] < 0)) {dev = Inf}
#  else if ((w[1] < 0) || (w[2] < 0) || (w[3] < 0) || (w[4] < 0)) {dev = Inf}
#  else {
#    datnew = dat %*% wnew
#    dev = scar(datnew,y,shape=c("in","in"))$deviance
#  }
#  return (dev)
# }
#
# index123 = optim(c(object$index[1:2,1],object$index[1:2,2]),fn)$par
# newindex = matrix(0,nrow=3,ncol=2)
# newindex[1:2,1] = index123[1:2]; newindex[1:2,2] = index123[3:4]
# newindex[3,1] = 1 - sum(abs(index123[1:2]))
# newindex[3,2] = 1 - sum(abs(index123[3:4]))
# newindex

Men's decathlon in 2012

Description

This dataset consists of men's decathlon results of at least 6500 points in 2012. Only those with athletes scored points in every event are included. Most men's decathlons are divided into a two-day competition. First day events include 100 metres, long jump, shot put, high jump and 400 metres. Second day events include 110 metres hurdles, discus throw, pole vault, javelin throw and 1500 metres.

Usage

data("decathlon_raw")

Format

A data frame containing 1075 observations on 15 variables.

Total

Total (decathlon) score.

Name

First name of the athlete.

Surname

Surname of the athlete.

X100m

Result of 100 metres, in seconds.

LJ

Result of long jump, in metres.

SP

Result of shot put, in metres.

HJ

Result of high jump, in metres.

X400m

Result of 400 metres, in seconds.

X110H

Result of 110 metres hurdles, in seconds.

DT

Result of discus throw, in metres.

PV

Result of polt vault, in metres.

JT

Result of javelin throw, in metres.

X1500m

Result of 1500 metres, in seconds.

First.day

Total points scored in the first day.

Second.day

Total points scored in the second day.

Source

Data compiled from the following sources:

www.decathlon2000.com/

www.iaaf.org/

www.ceskydesetiboj.wz.cz/statistika.html

See Also

decathlon

Plot additive index components of a "scair" object

Description

This function takes a fitted scair object produced by scair and plots the component functions of the additive indices that make it up, on the scale of the additive index predictors.

Usage

## S3 method for class 'scair'
plot(x, whichplot = 1:m, addp = TRUE, ...)

Arguments

Details

A plot is produced for each additive index component function in the fitted scair object. Each function is plotted in the range of the observed indices, with the identifiability condition satisfied (see Details of scar and scair).

Author(s)

Yining Chen and Richard Samworth

See Also

scair, predict.scair, plot.scar

Examples

## See examples for the function scair

Plot components of a "scar" object

Description

This function takes a fitted scar object produced by scar and plots the component functions that make it up, on the scale of the additive predictors.

Usage

## S3 method for class 'scar'
plot(x, whichplot = 1:d, addp = TRUE, ...)

Arguments

Details

A plot is produced for each additive component in the fitted scar object. Each function is plotted in the range of the observed covariates, with the identifiability condition satisfied (see Details of scar).

It can be invoked by calling "plot(x)" for an object x of the scar class, or directly by calling "plot.scar(x)" regardless of the class of the object.

Author(s)

Yining Chen and Richard Samworth

See Also

scar, predict.scar, plot.scair

Examples

## See examples for the function scar

Predict method for scair fits

Description

This function obtains predictions from a fitted scair object.

Usage

## S3 method for class 'scair'
predict(object, newdata, type = c("link", "response"), rule=1, ...)

Arguments

Value

A numeric vector of predictions.

Author(s)

Yining Chen and Richard Samworth

See Also

scair, plot.scair, predict.scar

Examples

## See examples for the function scair

Predict method for scar fits

Description

This function obtains predictions from a fitted scar object.

Usage

## S3 method for class 'scar'
predict(object, newdata, type = c("link", "response"), rule=1, ...)

Arguments

Value

A numeric vector of predictions.

Author(s)

Yining Chen and Richard Samworth

See Also

scar, plot.scar

Examples

## See examples for the function scar

Maximizing the likelihood of the generalised additive index model with shape constraints

Description

This function searches for a maximum likelihood estimator (mle) of the generalised additive index regression with shape constraints. A stochastic search strategy is used here.

Each index is a linear combination of some (or all) the covariates. Each additive component function of these index predictors is assumed to belong to one of the nine possible shape restrictions.

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

Usage

scair(x,y,shape=rep("l",1), family=gaussian(), weights=rep(1,length(y)), 
  epsilon=1e-8, delta=0.1, indexgen=c("unif", "norm"), iter = 200, 
  allnonneg = FALSE)

Arguments

Details

For i = 1,\ldots,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 index model can be written as

g(\mu) = f(x),

where x=(x_1,\ldots,x_d)^T, g is a known link function, A is an d \times m index matrix, and f is an additive function. Our task is to estimate both the index matrix and the additive function.

Assume the canonical link function is used here, then the maximum likelihood estimator of the generalised additive index model based on observations (X_1,Y_1), \ldots, (X_n,Y_n) is the function that maximises

\frac{1}{n} \sum_{i=1}^n w_i \{Y_i f(A^T X_i) - B(f(A^T X_i))\}

subject to the restrictions that for every j = 1,\ldots,m, 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.

For any given A, the optimization problem can solved using the active set algorithm implemented in scar. Therefore, this problem can be reduced to a finite-dimensional optimisation problem. Here we apply a simple stochastic search strategy is proposed, though other methods, such as downhill simplex, is also possible (and sometimes offers competitive performance). All the implementaton details can be found in Chen and Samworth (2016), where theoretical justification of our estimator (i.e. uniform consistency) is also given.

For the identifiability of additive index models, we refer to Yuan (2011).

Value

An object of class scair, with the following components:

Author(s)

Yining Chen and Richard Samworth

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.

Yuan, M. (2011). On the identifiability of additive index models. Statistica Sinica, 21, 1901-1911.

See Also

plot.scair, predict.scair, scar, decathlon

Examples

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

## Simulate the covariates and the responses
## covariates are drawn uniformly from [-1,1]^3
set.seed(10)
d = 3
n = 500
x = matrix(runif(n*d)*2-1,nrow=n,ncol=d) 
y = f(x) + rnorm(n,sd=0.5)

## Single index model so no delta is required here
shape=c("cvx")
object = scair(x,y,shape=shape, family=gaussian(),iter = 100)

## Estimated index matrix
object$index

## Root squared error for the estimated index
sqrt(sum((object$index - c(0.5,0.25,-0.25))^2))

## Plot the estimatied additive function for the single index
plot(object)

## Evaluate the estimated prediction 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))) 

## Here we can treat the obtained index matrix as a warm start and perform 
## further optimization (on the second and third entry of the index)
## using e.g. the default R optimisation routine.
fn<-function(w){
    dev = Inf
    if (abs(w[1])+abs(w[2])>1) return(dev)
    else {
      wnew = matrix(c(1-abs(w[1])-abs(w[2]),w[1],w[2]),ncol=1)
      dev = scar(x %*% wnew, y, shape = "cvx")$deviance
      return (dev)
    } 
}
index23 = optim(object$index[2:3],fn)$par
newindex = matrix(c(1-sum(abs(index23)),index23),ncol=1); newindex

## Root squared error for the new estimated index
sqrt(sum((newindex - c(0.5,0.25,-0.25))^2))

## A further example is provided in decathlon dataset

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

Details

For i = 1,\ldots,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(\mu) = f(x),

where x=(x_1,\ldots,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), \ldots, (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,\ldots,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,\ldots,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.

Value

An object of class scar, with the following components:

.

Note

We acknowledge that glm.fit from the R package stats is called to perform the method of iterated reweighted least squares (IRLS) in our routine. It is possible to speed up the implementation considerably by simply suppressing all the run-time checks there.

If all the component functions are linear, then it is prefered to call directly the function glm.

For the one-dimensional covariate, see the pool adjacent violators algorithm (PAVA) of Robertson, Wright and Dykstra (1998) and the support reduction method of Groeneboom, Jongbloed and Wellner (2008).

A different approach to tackle this problem is to use splines. See the R package scam. We stress here that our approach is free of tuning parameters while scam is not, which can be viewed as a major difference.

To estimate the generalised additive regression function without any shape restrictions, see Wood (2004) and Hastie and Tibshirani (1990). Their corresponding R implementations are mgcv and gam.

Author(s)

Yining Chen and Richard Samworth

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.

See Also

plot.scar, predict.scar, scair, scam, mgcv, gam, glm

Examples

## 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))) 

Internal shape-constrained additive regression functions

Description

Internal functions for the package scar

Details

These internal functions are not to be called by the user.

For a given real sequence a_1,..., a_n, the function aggsum returns the partial sums a_1, a_1+a_2, a_1+a_2+a_3,..., \sum_{i=1}^n a_i.