# cgam: Constrained Generalized Additive Model Fitting In cgam: Constrained Generalized Additive Model

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

 1 2 cgam(formula, nsim = 0, family = gaussian, cpar = 1.2, data = NULL, weights = NULL, sc_x = FALSE, sc_y = FALSE, pnt = TRUE, pen = 0, var.est = NULL) 

## 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 η 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 η and a predictor x in the following way: Assume that η is the systematic component and x is a predictor: incr(x): η is increasing in x. See incr for more details. s.incr(x): η is smoothly increasing in x. See s.incr for more details. decr(x): η is decreasing in x. See decr for more details. s.decr(x): η is smoothly decreasing in x. See s.decr for more details. conc(x): η is concave in x. See conc for more details. s.conc(x): η is smoothly concave in x. See s.conc for more details. conv(x): η is convex in x. See conv for more details. s.conv(x): η is smoothly convex in x. See s.conv for more details. incr.conc(x): η is increasing and concave in x. See incr.conc for more details. s.incr.conc(x): η is smoothly increasing and concave in x. See s.incr.conc for more details. decr.conc(x): η is decreasing and concave in x. See decr.conc for more details. s.decr.conc(x): η is smoothly decreasing and concave in x. See s.decr.conc for more details. incr.conv(x): η is increasing and convex in x. See incr.conv for more details. s.incr.conv(x): η is smoothly increasing and convex in x. See s.incr.conv for more details. decr.conv(x): η is decreasing and convex in x. See decr.conv for more details. s.decr.conv(x): η is smoothly decreasing and convex in x. See s.decr.conv for more details. s(x): η is smooth in x. See s for more details. tree(x): η has a tree-ordering in x. See tree for more details. umbrella(x): η 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 = 0. 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, 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 σ^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.

## Details

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

For the additive model, the systematic component is specified for each observation by η_i = f_1(x_{1i}) + … + f_L(x_{Li}) + z_i'β, where the functions f_l describe the relationships of the non-parametrically modelled predictors x_l, β 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 η is a n by 1 vector. The matrix form of the systematic component and the predictor is η = φ_1 + … + φ_L + Zβ, where φ_l is the individual component for the lth non-parametrically modelled predictor, l = 1, …, L, and Z is an n by p design matrix for the parametrically modelled covariate.

To model the component φ_l, smooth regression splines or non-smooth ordinal basis functions can be used. The constraints for the component φ_l are in C_l. In the first case, C_l = \{φ_l \in R^n: φ_l = v_l+B_lβ_l, where β_l ≥ 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 = \{φ \in R^n: A_lφ ≥ 0 and B_lφ = 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 + … + 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 η. muhat The fitted mean value, obtained by transforming the systematic component η 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 β. 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 β. A t-distribution is used as the approximate distribution. se.beta The standard errors for the estimation of the vector β. 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 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

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

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 # 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 Nov. 9, 2018, 1:04 a.m.