R语言含约束变量与二元变量的线性回归求解咨询

我正在运行一个包含3个自变量的线性回归模型：2个连续变量、1个代表节假日的二元变量，模型形式为：

log(Unit_Sales) ~ Intercept + log(Unit_Price) + log(Unit_Distribution) + Holiday（二元变量）

用lm()拟合后模型拟合度不错，但log(Unit_Distribution)的系数大于1，由于业务或统计需求，必须强制该系数小于1。尝试加入更多变量后该系数仍大于1，不想继续加变量，希望直接对系数加约束。

之前用optim()在排除二元变量时成功实现了约束，但加入Holiday后无法正常运行；用Rglpk_solve_LP时出现报错，推测是约束、rhs或上下界矩阵的写法有误。

请问该如何处理？有没有更合适的R包或函数？能否提供处理混合变量的Rglpk_solve_LP正确语法，或者其他解决方案？

解决方案 #

1. 用constrOptim快速实现带约束最小二乘 #

这是最直接的方法，不需要区分变量类型（连续/二元都能处理），核心是把最小二乘的目标函数定义为残差平方和，然后给目标系数加上限约束。

示例代码：

# 模拟你的数据集（替换成真实数据即可）
set.seed(123)
n <- 100
log_price <- rnorm(n)
log_dist <- rnorm(n)
holiday <- sample(c(0,1), n, replace = TRUE)
log_sales <- 0.5 + (-0.8)*log_price + 1.2*log_dist + 0.3*holiday + rnorm(n, 0, 0.2)
df <- data.frame(log_sales, log_price, log_dist, holiday)

# 定义最小二乘目标函数：计算残差平方和
ols_obj <- function(beta, X, y) {
  sum((y - X %*% beta)^2)
}

# 构造带截距项的设计矩阵
X <- model.matrix(~ log_price + log_dist + holiday, data = df)
y <- df$log_sales

# 用无约束回归的结果做初始值，加快收敛
init_beta <- coef(lm(log_sales ~ log_price + log_dist + holiday, data = df))

# 设置约束：log_dist的系数 ≤ 0.999（用0.999而非1避免边界数值问题）
# 约束矩阵：第3列对应log_dist的系数（需确认X的列顺序是否正确）
A <- matrix(0, nrow = 1, ncol = length(init_beta))
A[1, 3] <- 1
b <- 0.999

# 运行带约束的优化
constr_result <- constrOptim(theta = init_beta, 
                             f = ols_obj, 
                             grad = NULL,  # 自动计算数值梯度
                             ui = A, 
                             ci = b, 
                             X = X, 
                             y = y)

# 查看约束后的系数
constr_result$par

2. 用Rglpk_solve_QP做二次规划（你之前用错了函数） #

线性回归的最小二乘本质是二次优化问题，Rglpk_solve_LP是处理线性规划的，你应该用Rglpk_solve_QP。下面是正确的语法，同样支持混合变量：

library(Rglpk)

# 二次规划的核心矩阵：D = X^T X，d = X^T y
D <- t(X) %*% X
d <- t(X) %*% y

# 约束条件和之前一致：log_dist系数 ≤ 0.999
constraint_matrix <- A
constraint_rhs <- b

# 运行二次规划
qp_result <- Rglpk_solve_QP(Dmat = D, 
                            dvec = d, 
                            Amat = t(constraint_matrix),  # 注意这里要转置矩阵
                            bvec = constraint_rhs, 
                            meq = 0,  # 0表示所有约束都是不等式
                            lb = rep(-Inf, ncol(X)),  # 其他变量无下限
                            ub = rep(Inf, ncol(X)))  # 其他变量无上限

# 查看结果
qp_result$solution

3. 用glmnet做带约束的回归（适合兼顾正则化的场景） #

如果之后需要加正则化，glmnet的upper.limits参数可以直接给指定系数设上限，lambda设为0就是纯约束的最小二乘：

library(glmnet)

# 构造不含截距的矩阵（glmnet会自动加截距）
X_mat <- as.matrix(df[, c("log_price", "log_dist", "holiday")])
y_vec <- df$log_sales

# 设置各变量系数的上限：只有log_dist的上限是0.999，其他无限制
upper_limits <- c(Inf, 0.999, Inf)

# 拟合模型：alpha=0是岭回归，lambda=0关闭正则化，只应用约束
glmnet_result <- glmnet(X_mat, y_vec, 
                        alpha = 0, 
                        lambda = 0, 
                        upper.limits = upper_limits)

# 提取系数
coef(glmnet_result)

内容的提问来源于stack exchange，提问作者Roman Shuster