对系数有约束的线性回归
Linear regression with constraints on the coefficients
我正在尝试对这样的模型执行线性回归:
Y = aX1 + bX2 + c
所以,Y ~ X1 + X2
假设我有以下响应向量:
set.seed(1)
Y <- runif(100, -1.0, 1.0)
以及以下预测变量矩阵:
X1 <- runif(100, 0.4, 1.0)
X2 <- sample(rep(0:1,each=50))
X <- cbind(X1, X2)
我想对系数使用以下约束:
a + c >= 0
c >= 0
所以对 b 没有限制。
我知道 glmc 包可用于应用约束,但我无法确定如何为我的约束应用它。例如,我也知道可以使用 contr.sum 使所有系数之和为 0,但这不是我想要做的。 solve.QP() 似乎是另一种可能性,可以使用设置 meq=0
以便所有系数 >=0(同样,这不是我的目标)。
注意:解决方案必须能够处理响应向量 Y 中的 NA 值,例如:
Y <- runif(100, -1.0, 1.0)
Y[c(2,5,17,56,37,56,34,78)] <- NA
solve.QP
可以传递任意线性约束,因此它当然可以用于对约束 a+c >= 0
和 c >= 0
.
建模
首先,我们可以向 X
添加一列 1 以捕获截距项,然后我们可以使用 solve.QP
:
复制标准线性回归
X2 <- cbind(X, 1)
library(quadprog)
solve.QP(t(X2) %*% X2, t(Y) %*% X2, matrix(0, 3, 0), c())$solution
# [1] 0.08614041 0.21433372 -0.13267403
对于问题的示例数据,使用标准线性回归都不满足任何约束。
通过修改 Amat
和 bvec
参数,我们可以添加两个约束:
solve.QP(t(X2) %*% X2, t(Y) %*% X2, cbind(c(1, 0, 1), c(0, 0, 1)), c(0, 0))$solution
# [1] 0.0000000 0.1422207 0.0000000
根据这些约束,通过将 a 和 c 系数设置为都等于 0 来最小化平方残差。
您可以像 lm
函数一样处理 Y
或 X2
中的缺失值,方法是删除有问题的观察结果。作为预处理步骤,您可以执行类似以下操作:
has.missing <- rowSums(is.na(cbind(Y, X2))) > 0
Y <- Y[!has.missing]
X2 <- X2[!has.missing,]
我正在尝试对这样的模型执行线性回归:
Y = aX1 + bX2 + c
所以,Y ~ X1 + X2
假设我有以下响应向量:
set.seed(1)
Y <- runif(100, -1.0, 1.0)
以及以下预测变量矩阵:
X1 <- runif(100, 0.4, 1.0)
X2 <- sample(rep(0:1,each=50))
X <- cbind(X1, X2)
我想对系数使用以下约束:
a + c >= 0
c >= 0
所以对 b 没有限制。
我知道 glmc 包可用于应用约束,但我无法确定如何为我的约束应用它。例如,我也知道可以使用 contr.sum 使所有系数之和为 0,但这不是我想要做的。 solve.QP() 似乎是另一种可能性,可以使用设置 meq=0
以便所有系数 >=0(同样,这不是我的目标)。
注意:解决方案必须能够处理响应向量 Y 中的 NA 值,例如:
Y <- runif(100, -1.0, 1.0)
Y[c(2,5,17,56,37,56,34,78)] <- NA
solve.QP
可以传递任意线性约束,因此它当然可以用于对约束 a+c >= 0
和 c >= 0
.
首先,我们可以向 X
添加一列 1 以捕获截距项,然后我们可以使用 solve.QP
:
X2 <- cbind(X, 1)
library(quadprog)
solve.QP(t(X2) %*% X2, t(Y) %*% X2, matrix(0, 3, 0), c())$solution
# [1] 0.08614041 0.21433372 -0.13267403
对于问题的示例数据,使用标准线性回归都不满足任何约束。
通过修改 Amat
和 bvec
参数,我们可以添加两个约束:
solve.QP(t(X2) %*% X2, t(Y) %*% X2, cbind(c(1, 0, 1), c(0, 0, 1)), c(0, 0))$solution
# [1] 0.0000000 0.1422207 0.0000000
根据这些约束,通过将 a 和 c 系数设置为都等于 0 来最小化平方残差。
您可以像 lm
函数一样处理 Y
或 X2
中的缺失值,方法是删除有问题的观察结果。作为预处理步骤,您可以执行类似以下操作:
has.missing <- rowSums(is.na(cbind(Y, X2))) > 0
Y <- Y[!has.missing]
X2 <- X2[!has.missing,]