R: Shape-Restricted Regression

shapereg {coneproj}

R Documentation

Shape-Restricted Regression

Description

The regression model y_i = f(t_i) + x_i'\beta + \varepsilon_i, i = 1,\ldots,n is considered, where the only assumptions about f concern its shape. The vector expression for the model is y = \theta + X\beta + \varepsilon. X represents a parametrically modelled covariate, which could be a categorical covariate or a linear term. The shapereg function allows eight shapes: increasing, decreasing, convex, concave, increasing-convex, increasing-concave, decreasing-convex, and decreasing-concave. This routine employs a single cone projection to find \theta and \beta simultaneously.

Usage

shapereg(formula, data = NULL, weights = NULL, test = FALSE, nloop = 1e+4)

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`. A predictor can be a non-parametrically modelled variable with a shape restriction or a parametrically modelled unconstrained covariate. In terms of a non-parametrically modelled predictor, the user is supposed to indicate the relationship between `E(y)` and a predictor `t` in the following way: incr(t): `E(y)` is increasing in `t`. See `incr` for more details. decr(t): `E(y)` is decreasing in `t`. See `decr` for more details. conc(t): `E(y)` is concave in `t`. See `conc` for more details. conv(t): `E(y)` is convex in `t`. See `conv` for more details. incr.conc(t): `E(y)` is increasing and concave in `t`. See `incr.conc` for more details. decr.conc(t): `E(y)` is decreasing and concave in `t`. See `decr.conc` for more details. incr.conv(t): `E(y)` is increasing and convex in `t`. See `incr.conv` for more details. decr.conv(t): `E(y)` is decreasing and convex in `t`. See `decr.conv` 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.
`test`	The test parameter given by the user.
`nloop`	The number of simulations used to get the p-value for the `E_{01}` test. The default is 1e+4.

Details

This routine constrains \theta in the equation y = \theta + X\beta + \varepsilon by a shape parameter.

The constraint cone C has the form \{\phi: \phi = v + \sum b_i\delta_i, i = 1,\ldots,m, b_1,\ldots, b_m \ge 0 \}, v is in V. The column vectors of X are in V, i.e., the linear space contained in the constraint cone.

The hypothesis test H_0: \phi is in V versus H_1: \phi is in C is an exact one-sided test, and the test statistic is E_{01} = (SSE_0 - SSE_1)/(SSE_0), which has a mixture-of-betas distribution when H_0 is true and \varepsilon is a vector following a standard multivariate normal distribution with mean 0. The mixing parameters are found through simulations. The number of simulations used to obtain the mixing distribution parameters for the test is 10,000. Such simulations usually take some time. For the "feet" data set used as an example in this section, whose sample size is 39, the time to get a p-value is roughly between 4 seconds.

This routine calls coneB for the cone projection part.

Value

`coefs`	The estimated coefficients for `X`, i.e., the estimation for the vector `\beta`. Note that even if the user does not provide a constant vector in `X`, the coefficient for the intercept will be returned.
`constr.fit`	The shape-restricted fit over the constraint cone `C` of the form `\{\phi: \phi = v + \sum b_i\delta_i, i = 1,\ldots,m, b_1,\ldots, b_m \ge 0 \}`, `v` is in `V`.
`linear.fit`	The least-squares regression of `y` on `V`, i.e., the linear space contained in the constraint cone. If shape is 3 or shape is 4, `V` is spanned by `X` and `t`. Otherwise, it is spanned by `X`. `X` must be full column rank, and the matrix formed by combining `X` and `t` must also be full column rank.
`se.beta`	The standard errors for the estimation of the vector `\beta`. The degree of freedom is returned by coneB and is multiplied by 1.5. Note that even if the user does not provide a constant vector in `X`, the standard error for the intercept will be returned.
`pval`	The p-value for the hypothesis test `H_0: \phi` is in `V` versus `H_1: \phi` is in `C`. `C` is the constraint cone of the form `\{\phi: \phi = v + \sum b_i\delta_i, i = 1,\ldots,m, b_1,\ldots, b_m \ge 0 \}`, `v` is in `V`, and `V` is the linear space contained in the constraint cone. If test == TRUE, a p-value is returned. Otherwise, the test is skipped and no p-value is returned.
`pvals.beta`	The approximate p-values for the estimation of the vector `\beta`. A t-distribution is used as the approximate distribution. Note that even if the user does not provide a constant vector in `X`, the approximate p-value for the intercept will be returned.
`test`	The test parameter given by the user.
`SSE0`	The sum of squared residuals for the linear part.
`SSE1`	The sum of squared residuals for the full model.
`shape`	A number showing the shape constraint given by the user in a shapereg fit.
`tms`	The terms objects extracted by the generic function terms from a shapereg fit.
`zid`	A vector keeping track of the position 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.
`tnm`	The name of the shape-restricted predictor.
`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 factor or a linear term.
`is_fac`	A logical scalar showing if or not a variable is a factor.
`xmat`	A matrix whose columns represent the parametrically modelled covariate.
`call`	The matched call.

Author(s)

Mary C. Meyer and Xiyue Liao

References

Raubertas, R. F., C.-I. C. Lee, and E. V. Nordheim (1986) Hypothesis tests for normals means constrained by linear inequalities. Communications in Statistics - Theory and Methods 15 (9), 2809–2833.

Robertson, T., F. Wright, and R. Dykstra (1988) Order Restricted Statistical Inference New York: John Wiley and Sons.

Fraser, D. A. S. and H. Massam (1989) A mixed primal-dual bases algorithm for regression under inequality constraints application to concave regression. Scandinavian Journal of Statistics 16, 65–74.

Meyer, M. C. (2003) A test for linear vs convex regression function using shape-restricted regression. Biometrika 90(1), 223–232.

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

Meyer, M.C.(2013a) Semiparametric additive constrained regression. Journal of Nonparametric Statistics 25(3), 715–743.

Liao, X. and M. C. Meyer (2014) coneproj: An R package for the primal or dual cone projections with routines for constrained regression. Journal of Statistical Software 61(12), 1–22.

Examples

# load the feet data set
    data(feet)

# extract the continuous and constrained predictor
    l <- feet$length

# extract the continuous response
    w <- feet$width

# extract the categorical covariate: sex
    s <- feet$sex

# make an increasing fit with test set as FALSE
    ans <- shapereg(w ~ incr(l) + factor(s))

# check the summary table 
    summary(ans)

# make an increasing fit with test set as TRUE
    ans <- shapereg(w ~ incr(l) + factor(s), test = TRUE, nloop = 1e+3)

# check the summary table 
    summary(ans)

# make a plot comparing the unconstrained fit and the constrained fit
    par(mar = c(4, 4, 1, 1))
    ord <- order(l)
    plot(sort(l), w[ord], type = "n", xlab = "foot length (cm)", ylab = "foot width (cm)")
    title("Shapereg Example Plot")

# sort l according to sex
    ord1 <- order(l[s == "G"])
    ord2 <- order(l[s == "B"])

# make the scatterplot of l vs w for boys and girls
    points(sort(l[s == "G"]), w[s == "G"][ord1], pch = 21, col = 1)
    points(sort(l[s == "B"]), w[s == "B"][ord2], pch = 24, col = 2)

# make an unconstrained fit to boys and girls
    fit <- lm(w ~ l + factor(s))

# plot the unconstrained fit 
    lines(sort(l), (coef(fit)[1] + coef(fit)[2] * l + coef(fit)[3])[ord], lty = 2)
    lines(sort(l), (coef(fit)[1] + coef(fit)[2] * l)[ord], lty = 2, col = 2)
    legend(21.5, 9.8, c("boy","girl"), pch = c(24, 21), col = c(2, 1)) 

# plot the constrained fit
    lines(sort(l), (ans$constr.fit - ans$linear.fit + coef(ans)[1])[ord], col = 1)
    lines(sort(l), (ans$constr.fit - ans$linear.fit + coef(ans)[1] + coef(ans)[2])[ord], col = 2)

[Package coneproj version 1.19 Index]