constreg function

Constrained Parametric Regression

The least-squares regression model $y = X\beta + \varepsilon$ is considered, where the object is to find $\beta$ to minimize $||y - X\beta||^2$, subject to $A\beta \ge 0$.

constreg(y, xmat, amat, w = NULL, test = FALSE, nloop = 1e+4)

Arguments

• y: A vector of length $n$.
• xmat: A full column-rank design matrix. The column number of xmat must equal the length of $\beta$.
• amat: A constraint matrix. The rows of amat must be irreducible. The column number of amat must equal the length of $\beta$.
• w: An optional nonnegative vector of weights of length $n$. If w is not given, all weights are taken to equal 1. Otherwise, the minimization of $(y - X\beta)'w(y - X\beta)$ over $C$ is returned. The default is w = NULL.
• test: A logical scalar. If test == TRUE, then the p-value for the test $H_0:\beta$ is in $V$ versus $H_1:\beta$ is in $C$ is returned. $C$ is the constraint cone of the form $\{\beta: A\beta \ge 0\}$, and $V$ is the null space of $A$. The default is test = FALSE.
• nloop: The number of simulations used to get the p-value for the $E_{01}$ test. The default is 1e+4.

Returns

• constr.fit: The constrained fit of $y$ given that $\beta$ is in the cone $C$ of the form $\{\beta: A\beta \ge 0 \}$.

• unconstr.fit: The unconstrainted fit, i.e., the least-squares regression of $y$ on the space spanned by $X$.

• pval: The p-value for the hypothesis test $H_0:\beta$ is in $V$ versus $H_1:\beta$ is in $C$. The constraint cone $C$ has the form $\{\beta: A\beta \ge 0 \}$ and $V$ is the null space of $A$. If test == TRUE, a p-value is returned. Otherwise, the test is skipped and no p-value is returned.

• coefs: The estimated constrained parameters, i.e., the estimation of the vector $\beta$.

Details

The hypothesis test $H_0:\beta$ is in $V$ versus $H_1:\beta$ 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 "FEV" data set used as an example in this section, whose sample size is 654, the time to get a p-value is roughly 6 seconds.

The constreg function calls coneA for the cone projection part.

References

Brunk, H. D. (1958) On the estimation of parameters restricted by inequalities. The Annals of Mathematical Statistics 29 (2), 437--454.

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.

Meyer, M. C. and J. C. Wang (2012) Improved power of one-sided tests. Statistics and Probability Letters 82, 1619--1622.

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.

Author(s)

Mary C. Meyer and Xiyue Liao

See Also

coneA

Examples

# load the FEV data set
    data(FEV)

# extract the variables
    y <- FEV$FEV
    age <- FEV$age
    height <- FEV$height
    sex <- FEV$sex
    smoke <- FEV$smoke

# scale age and height
    scale_age <- (age - min(age)) / (max(age) - min(age))
    scale_height <- (height - min(height)) / (max(height) - min(height))

# make xmat
    xmat <- cbind(1, scale_age, scale_height, scale_age * scale_height, sex, smoke)

# make the constraint matrix 
    amat <- matrix(0, 4, 6)
    amat[1, 2] <- 1; amat[2, 2] <- 1; amat[2, 4] <- 1 
    amat[3, 3] <- 1; amat[4, 3] <- 1; amat[4, 4] <- 1

# call constreg to get constrained coefficient estimates
    ans1 <- constreg(y, xmat, amat)
    bhat1 <- coef(ans1)

# call lm to get unconstrained coefficient estimates
    ans2 <- lm(y ~ xmat[,-1])
    bhat2 <- coef(ans2)

# create a 3D plot to show the constrained fit and the unconstrained fit 
    n <- 25
    xgrid <- seq(0, 1, by = 1/n)
    ygrid <- seq(0, 1, by = 1/n)
    x1 <- rep(xgrid, each = length(ygrid))
    x2 <- rep(ygrid, length(xgrid))
    xinterp <- cbind(x1, x2)
    xmatp <- cbind(1, xinterp, x1 * x2, 0, 0)
    
    thint1 <- crossprod(t(xmatp), bhat1)
    A1 <- matrix(thint1, length(xgrid), length(ygrid), byrow = TRUE) 
    thint2 <- crossprod(t(xmatp), bhat2)
    A2 <- matrix(thint2, length(xgrid), length(ygrid), byrow = TRUE) 

    par(mfrow = c(1, 2))
    par(mar = c(4, 1, 1, 1))
    persp(xgrid, ygrid, A1, xlab = "age", ylab = "height", 
    zlab = "FEV", theta = -30)
    title("Constrained Fit")

    par(mar = c(4, 1, 1, 1))
    persp(xgrid, ygrid, A2, xlab = "age", ylab = "height", 
    zlab = "FEV", theta = -30)
    title("Unconstrained Fit")

Note

In the 3D plot of the "FEV" example, it is shown that the unconstrained fit increases as "age" increases when "height" is large, but decreases as "age" increases when "height" is small. This does not make sense, since "FEV" should not decrease with respect to "age" given any value of "height". The constrained fit avoids this situation by keeping the fit of "FEV" non-decreasing with respect to "age".