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, cic = FALSE, nsim = 100, family = gaussian, cpar = 1.5, 
	data = NULL, weights = NULL, sc_x = FALSE, sc_y = FALSE, pnt = TRUE, 
	pen = 0, var.est = NULL, gcv = FALSE, pvf = TRUE)

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.
edf: The constrained effective degrees of freedom.
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.
y: The response variable.
xmat_add: 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.
prior.w: User-defined weights.
wt: The weights in the final iteration of the iteratively re-weighted cone projections.
wt.iter: Logical flag 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.
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 names of the shape-restricted predictors and the smooth predictors with no shape constraint in xmat.
ynm: The name of the response variable.
znms: A vector storing the names 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
vh: The estimated monotonic variance function.
kts.var: The knots used in monotonic variance function estimation.
call: The matched call.

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.

cic

Logical flag indicating if or not simulations are used to get the cic value. The default is cic = FALSE

nsim

The number of simulations used to get the cic parameter. The default is nsim = 100.

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 four families used in csvy: Gaussian, binomial, poisson, and Gamma. Note that if family = Ord is specified, a proportional odds regression model with shape constraints is fitted. This is under development.

cpar

A multiplier to estimate the model variance, which is defined as \(\sigma^2 = SSR / (n - cpar * edf)\). SSR is the sum of squared residuals for the full model and edf is the effective degrees of freedom. The default is cpar = 1.2. The user-defined value must be between 1 and 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

Logical flag indicating if or not continuous predictors are normalized. The default is sc_x = FALSE.

sc_y

Logical flag indicating if or not the response variable is normalized. The default is sc_y = FALSE.

pen

User-defined penalty parameter. It must be non-negative. It will only be used in a warped-plane spline fit or a triangle spline fit. The default is pen = 0.

pnt

Logical flag indicating if or not penalized constrained regression splines are used. It will only be used in a warped-plane spline fit or a triangle spline fit. The default is pnt = TRUE.

var.est

To do a monotonic variance function estimation, the user can set var.est = s.incr(x) or var.est = s.decr(x). See s.incr and s.decr for more details. The default is var.est = NULL.

gcv

Logical flag indicating if or not gcv is used to choose a penalty term in warped-plane surface fit. The default is gcv = FALSE.

pvf

Logical flag indicating if or not simulations are used to find the p-value of the test of linear vs double monotone in warped plane surface fit.

Author

Mary C. Meyer and Xiyue Liao

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

Liao, X. and Meyer, M. C. (2019) cgam: An R Package for the Constrained Generalized Additive Model. Journal of Statistical Software 89(5), 1--24.

Meyer, M. C. (2018) A Framework for Estimation and Inference in Generalized Additive Models with Shape and Order Restrictions. Statistical Science 33(4), 595--614.

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.
if (FALSE) {
  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)
  summary(fit.cgam.s)
  anova(fit.cgam.s)
  
  # 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")

# Example 4. monotonic variance estimation
  n <- 400
  x <- runif(n)
  sig <- .1 + exp(15*x-8)/(1+exp(15*x-8))
  e <- rnorm(n)
  mu <- 10*x^2
  y <- mu + sig*e

  fit <- cgam(y ~ s.incr.conv(x), var.est = s.incr(x))
  est.var <- fit$vh
  muhat <- fit$muhat

  par(mfrow = c(1, 2))
  plot(x, y)
  points(sort(x), muhat[order(x)], type = "l", lwd = 2, col = 2)
  lines(sort(x), (mu)[order(x)], col = 4)

  plot(sort(x), est.var[order(x)], col=2, lwd=2, type="l", 
  lty=2, ylab="Variance", ylim=c(0, max(c(est.var, sig^2))))
  points(sort(x), (sig^2)[order(x)], col=1, lwd=2, type="l")

# Example 5. monotonic variance estimation with the lidar data set in SemiPar
  library(SemiPar)
  data(lidar)

  fit <- cgam(logratio ~ s.decr(range), var.est=s.incr(range), data=lidar)
  muhat <- fit$muhat
  est.var <- fit$vh
  
  logratio <- lidar$logratio
  range <- lidar$range
  pfit <- predict(fit, newData=data.frame(range=range), interval="confidence", level=0.95)
  upp <- pfit$upper
  low <- pfit$lower
  
  par(mfrow = c(1, 2))
  plot(range, logratio)
  points(sort(range), muhat[order(range)], type = "l", lwd = 2, col = 2)
  lines(sort(range), upp[order(range)], type = "l", lwd = 2, col = 4)
  lines(sort(range), low[order(range)], type = "l", lwd = 2, col = 4)
  title("Smoothly Decreasing Fit with a Point-Wise Confidence Interval", cex.main=0.5)
  
  plot(range, est.var, col=2, lwd=2, type="l",lty=2, ylab="variance")
  title("Smoothly Increasing Variance", cex.main=0.5)

Run the code above in your browser using DataLab