cgam: Constrained Generalized Additive Model Fitting

In cgam: Constrained Generalized Additive Model

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)

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.

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 lth 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.

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.

Author(s)

Mary C. Meyer and Xiyue Liao

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

# 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.
## Not run: 
  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) 

## End(Not run)

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

cgam documentation built on Aug. 10, 2023, 5:11 p.m.