Learn R Programming

cgam (version 1.6)

cgam: Constrained Generalized Additive Model Fitting

Description

The partial linear generalized additive model is fitted using the method of maximum likelihood, where shape or order restrictions can be imposed on the non-parametrically modelled predictors with optional smoothing, and no restrictions are imposed on the optional parametrically modelled covariate.

Usage

cgam(formula, nsim = 1e+2, family = gaussian(), cpar = 1.2, data = NULL, weights = NULL,
sc_x = FALSE, sc_y = FALSE)

Arguments

formula

A formula object which gives a symbolic description of the model to be fitted. It has the form "response ~ predictor". The response is a vector of length \(n\). The specification of the model can be one of the three exponential families: gaussian, binomial and poisson. The systematic component \(\eta\) is \(E(y)\), the log odds of \(y = 1\), and the logarithm of \(E(y)\) respectively. A predictor can be a non-parametrically modelled variable with or without a shape or order restriction, or a parametrically modelled unconstrained covariate. In terms of a non-parametrically modelled predictor, the user is supposed to indicate the relationship between the systematic component \(\eta\) and a predictor \(x\) in the following way:

Assume that \(\eta\) is the systematic component and \(x\) is a predictor:

  • incr(x): \(\eta\) is increasing in \(x\). See incr for more details.

  • s.incr(x): \(\eta\) is smoothly increasing in \(x\). See s.incr for more details.

  • decr(x): \(\eta\) is decreasing in \(x\). See decr for more details.

  • s.decr(x): \(\eta\) is smoothly decreasing in \(x\). See s.decr for more details.

  • conc(x): \(\eta\) is concave in \(x\). See conc for more details.

  • s.conc(x): \(\eta\) is smoothly concave in \(x\). See s.conc for more details.

  • conv(x): \(\eta\) is convex in \(x\). See conv for more details.

  • s.conv(x): \(\eta\) is smoothly convex in \(x\). See s.conv for more details.

  • incr.conc(x): \(\eta\) is increasing and concave in \(x\). See incr.conc for more details.

  • s.incr.conc(x): \(\eta\) is smoothly increasing and concave in \(x\). See s.incr.conc for more details.

  • decr.conc(x): \(\eta\) is decreasing and concave in \(x\). See decr.conc for more details.

  • s.decr.conc(x): \(\eta\) is smoothly decreasing and concave in \(x\). See s.decr.conc for more details.

  • incr.conv(x): \(\eta\) is increasing and convex in \(x\). See incr.conv for more details.

  • s.incr.conv(x): \(\eta\) is smoothly increasing and convex in \(x\). See s.incr.conv for more details.

  • decr.conv(x): \(\eta\) is decreasing and convex in \(x\). See decr.conv for more details.

  • s.decr.conv(x): \(\eta\) is smoothly decreasing and convex in \(x\). See s.decr.conv for more details.

  • s(x): \(\eta\) is smooth in \(x\). See s for more details.

  • tree(x): \(\eta\) has a tree-ordering in \(x\). See tree for more details.

  • umbrella(x): \(\eta\) has an umbrella-ordering in \(x\). See umbrella for more details.

nsim

The number of simulations used to get the cic parameter. The default is nsim = 1e+2.

family

A parameter indicating the error distribution and link function to be used in the model. It can be a character string naming a family function or the result of a call to a family function. This is borrowed from the glm routine in the stats package. There are three families used in cgam: gaussian, binomial and poisson.

cpar

A multiplier to estimate the model variance, which is defined as \(\sigma^2 = SSR / (n - d_0 - cpar * edf)\). SSR is the sum of squared residuals for the full model, \(d_0\) is the dimension of the linear space in the cone, and edf is the effective degrees of freedom. The default is cpar = 1.2. See Meyer, M. C. and M. Woodroofe (2000) for more details.

data

An optional data frame, list or environment containing the variables in the model. The default is data = NULL.

weights

An optional non-negative vector of "replicate weights" which has the same length as the response vector. If weights are not given, all weights are taken to equal 1. The default is weights = NULL.

sc_x

A logical scalar showing if or not non-parametrically modelled predictors will be scaled. The default is sc_x = FALSE.

sc_y

A logical scalar showing if or not the response will be scaled. The default is sc_y = FALSE.

Value

etahat

The fitted systematic component \(\eta\).

muhat

The fitted mean value, obtained by transforming the systematic component \(\eta\) by the inverse of the link function.

vcoefs

The estimated coefficients for the basis spanning the null space of the constraint set.

xcoefs

The estimated coefficients for the edges corresponding to the smooth predictors with no shape constraint and shape-restricted predictors.

zcoefs

The estimated coefficients for the parametrically modelled covariate, i.e., the estimation for the vector \(\beta\).

ucoefs

The estimated coefficients for the edges corresponding to the predictors with an umbrella-ordering constraint.

tcoefs

The estimated coefficients for the edges corresponding to the predictors with a tree-ordering constraint.

coefs

The estimated coefficients for the basis spanning the null space of the constraint set and edges corresponding to the shape-restricted and order-restricted predictors.

cic

The cone information criterion proposed in Meyer(2013a). It uses the "null expected degrees of freedom" as a measure of the complexity of the model. See Meyer(2013a) for further details of cic.

d0

The dimension of the linear space contained in the cone generated by all constraint conditions.

edf0

The estimated "null expected degrees of freedom". It is a measure of the complexity of the model. See Meyer (2013a) and Meyer (2013b) for further details.

etacomps

The fitted systematic component value for non-parametrically modelled predictors. It is a matrix of which each row is the fitted systematic component value for a non-parametrically modelled predictor. If there are more than one such predictors, the order of the rows is the same as the order that the user defines such predictors in the formula argument of cgam.

xmat

A matrix whose columns represent the shape-restricted predictors and smooth predictors with no shape constraint.

zmat

A matrix whose columns represent the basis for the parametrically modelled covariate. The user can choose to include a constant vector in it or not. It must have full column rank.

ztb

A list keeping track of the order of the parametrically modelled covariate.

tr

A matrix whose columns represent the predictors with a tree-ordering constraint.

umb

A matrix whose columns represent the predictors with an umbrella-ordering constraint.

tree.delta

A matrix whose rows are the edges corresponding to the predictors with a tree-ordering constraint.

umbrella.delta

A matrix whose rows are the edges corresponding to the predictors with an umbrella-ordering constraint.

bigmat

A matrix whose rows are the basis spanning the null space of the constraint set and the edges corresponding to the shape-restricted and order-restricted predictors.

shapes

A vector including the shape and partial-ordering constraints in a cgam fit.

shapesx

A vector including the shape constraints in a cgam fit.

wt

The weights in the final iteration of the iteratively re-weighted cone projections.

wt.iter

A logical scalar indicating if or not iteratively re-weighted cone projections may be used. If the response is gaussian, then wt.iter = FALSE; if the response is binomial or poisson, then wt.iter = TRUE.

family

The family parameter defined in a cgam formula.

SSE0

The sum of squared residuals for the linear part.

SSE1

The sum of squared residuals for the full model.

pvals.beta

The approximate p-values for the estimation of the vector \(\beta\). A t-distribution is used as the approximate distribution.

se.beta

The standard errors for the estimation of the vector \(\beta\).

null_df

The degree of freedom for the null model of a cgam fit, i.e., the model only containing a constant vector.

df

The degree of freedom for the null space of a cgam fit.

resid_df_obs

The observed degree of freedom for the residuals of a cgam fit.

null_deviance

The deviance for the null model of a cgam fit, i.e., the model only containing a constant vector.

deviance

The residual deviance of a cgam fit.

tms

The terms objects extracted by the generic function terms from a cgam fit. See the official help page (http://stat.ethz.ch/R-manual/R-patched/library/stats/html/terms.html) of the terms function for more details.

capm

The number of edges corresponding to the shape-restricted predictors.

capms

The number of edges corresponding to the smooth predictors with no shape constraint.

capk

The number of non-constant columns of zmat.

capt

The number of edges corresponding to the tree-ordering predictors.

capu

The number of edges corresponding to the umbrella-ordering predictors.

xid1

A vector keeping track of the beginning position of the set of edges in bigmat for each shape-restricted predictor and smooth predictor with no shape constraint in xmat.

xid2

A vector keeping track of the end position of the set of edges in bigmat for each shape-restricted predictor and smooth predictor with no shape constraint in xmat.

tid1

A vector keeping track of the beginning position of the set of edges in bigmat for each tree-ordering factor in tr.

tid2

A vector keeping track of the end position of the set of edges in bigmat for each tree-ordering factor in tr.

uid1

A vector keeping track of the beginning position of the set of edges in bigmat for each umbrella-ordering factor in umb.

uid2

A vector keeping track of the end position of the set of edges in bigmat for each umbrella-ordering factor in umb.

zid

A vector keeping track of the positions of the parametrically modelled covariate.

vals

A vector storing the levels of each variable used as a factor.

zid1

A vector keeping track of the beginning position of the levels of each variable used as a factor.

zid2

A vector keeping track of the end position of the levels of each variable used as a factor.

nsim

The number of simulations used to get the cic parameter.

xnms

A vector storing the name of the shape-restricted predictor and the smooth predictor with no shape constraint in xmat.

ynm

The name of the response variable.

znms

A vector storing the name of the parametrically modelled covariate.

is_param

A logical scalar showing if or not a variable is a parametrically modelled covariate, which could be a linear term or a factor.

is_fac

A logical scalar showing if or not a variable is a factor.

knots

A list storing the knots used for each shape-restricted predictor and smooth predictor with no shape constraint. For a smooth, constrained and a smooth, unconstrainted predictor, knots is a vector of more than \(1\) elements, and for a shape-restricted predictor without smoothing, knots = \(0\).

numknots

A vector storing the number of knots for each shape-restricted predictor and smooth predictor with no shape constraint. For a smooth, constrained and a smooth, unconstrainted predictor, numknots > \(1\), and for a shape-restricted predictor without smoothing, numknots = \(0\).

sps

A character vector storing the space parameter to create knots for each shape-restricted predictor.

ms

The centering terms used to make edges for shape-restricted predictors.

cpar

The cpar argument in the cgam formula

pl

A vector of the placebo values for the non-parametrically modelled covariates with a tree-ordering constraint.

call

The matched call.

Details

We consider generalized partial linear models with independent observations from an exponential family of the form \(p(y_i;\theta,\tau) = exp[\{y_i\theta_i - b(\theta_i)\}\tau - c(y_i, \tau)], i = 1,\ldots,n\), where the specifications of the functions \(b\) and \(c\) determine the sub-family of models. The mean vector \(\mu = E(y)\) has values \(\mu_i = b'(\theta_i)\), and is related to a design matrix of predictor variables through a monotonically increasing link function \(g(\mu_i) = \eta_i, i = 1,\ldots,n\), where \(\eta\) is the systematic component and describes the relationship with the predictors. The relationship between \(\eta\) and \(\theta\) is determined by the link function \(b\).

For the additive model, the systematic component is specified for each observation by \(\eta_i = f_1(x_{1i}) + \ldots + f_L(x_{Li}) + z_i'\beta\), where the functions \(f_l\) describe the relationships of the non-parametrically modelled predictors \(x_l\), \(\beta\) is a parameter vector, and \(z_i\) contains the values of variables to be modelled parametrically. The non-parametric components are modelled with shape or order assumptions with optional smoothing, and the solution is obtained through an iteratively re-weighted cone projection, with no back-fitting of individual components.

Suppose that \(\eta\) is a \(n\) by \(1\) vector. The matrix form of the systematic component and the predictor is \(\eta = \phi_1 + \ldots + \phi_L + Z\beta\), where \(\phi_l\) is the individual component for the \(l\)th non-parametrically modelled predictor, \(l = 1, \ldots, L\), and \(Z\) is an \(n\) by \(p\) design matrix for the parametrically modelled covariate.

To model the component \(\phi_l\), smooth regression splines or non-smooth ordinal basis functions can be used. The constraints for the component \(\phi_l\) are in \(C_l\). In the first case, \(C_l\) = \(\{\phi_l \in R^n: \phi_l = v_l+B_l\beta_l\), where \(\beta_l \ge 0\) and \(v_l\in V_l \}\), where \(B_l\) has regression splines as columns and \(V_l\) is the linear space contained in \(C_l\), and in the second case, \(C_l\) = \(\{\phi \in R^n: A_l\phi \ge 0\) and \(B_l\phi = 0\}\), for inequality constraint matrix \(A_l\) and equality constraint matrix \(B_l\).

The set \(C_l\) is a convex cone and the set \(C = C_1 + \ldots + C_p + Z\) is also a convex cone with a finite set of edges, where the edges are the generators of \(C\), and \(Z\) is the column space of the design matrix \(Z\) for the parametrically modelled covariate.

An iteratively re-weighted cone projection algorithm is used to fit the generalized regression model over the cone \(C\).

See references cited in this section and the official manual (https://cran.r-project.org/package=coneproj) for the R package coneproj for more details.

References

Meyer, M. C. (2013a) Semi-parametric additive constrained regression. Journal of Nonparametric Statistics 25(3), 715

Meyer, M. C. (2013b) A simple new algorithm for quadratic programming with applications in statistics. Communications in Statistics 42(5), 1126--1139.

Meyer, M. C. and M. Woodroofe (2000) On the degrees of freedom in shape-restricted regression. Annals of Statistics 28, 1083--1104.

Meyer, M. C. (2008) Inference using shape-restricted regression splines. Annals of Applied Statistics 2(3), 1013--1033.

Mammen, E. and K. Yu (2007) Additive isotonic regression. IMS Lecture Notes-Monograph Series Asymptotics: Particles, Process, and Inverse Problems 55, 179--195.

Huang, J. (2002) A note on estimating a partly linear model under monotonicity constraints. Journal of Statistical Planning and Inference 107, 343--351.

Cheng, G.(2009) Semiparametric additive isotonic regression. Journal of Statistical Planning and Inference 139, 1980--1991.

Bacchetti, P. (1989) Additive isotonic models. Journal of the American Statistical Association 84(405), 289--294.

Examples

Run this code
# Example 1.
  data(cubic)
  # extract x
  x <- cubic$x

  # extract y
  y <- cubic$y

  # regress y on x with no restriction with lm()
  fit.lm <- lm(y ~ x + I(x^2) + I(x^3))

  # regress y on x under the restriction: "increasing and convex"
  fit.cgam <- cgam(y ~ incr.conv(x))

  # make a plot to compare the two fits
  par(mar = c(4, 4, 1, 1))
  plot(x, y, cex = .7, xlab = "x", ylab = "y")
  lines(x, fit.cgam$muhat, col = 2, lty = 2)
  lines(x, fitted(fit.lm), col = 1, lty = 1)
  legend("topleft", bty = "n", c("constrained cgam fit", "unconstrained lm fit"), 
  lty = c(2, 1), col = c(2, 1))

# Example 2.

  library(gam)
  data(kyphosis)
  
  # regress Kyphosis on Age, Number, and Start under the restrictions:
  # "concave", "increasing and concave", and "decreasing and concave" 
  fit <- cgam(Kyphosis ~ conc(Age) + incr.conc(Number) + decr.conc(Start), 
  family = binomial(), data = kyphosis) 


# Example 3.
  library(MASS)
  data(Rubber)
  
  # regress loss on hard and tens under the restrictions:
  # "decreasing" and "decreasing"
  fit.cgam <- cgam(loss ~ decr(hard) + decr(tens), data = Rubber)
  # "smooth and decreasing" and "smooth and decreasing"
  fit.cgam.s <- cgam(loss ~ s.decr(hard) + s.decr(tens), data = Rubber)

  # make a 3D plot based on fit.cgam and fit.cgam.s
  plotpersp(fit.cgam, th = 120, main = "3D Plot of a Cgam Fit")
  plotpersp(fit.cgam.s, tens, hard, data = Rubber, th = 120, main = "3D Plot of a Smooth Cgam Fit")

Run the code above in your browser using DataLab