Линейная регрессия с ограничениями на коэффициенты

Question

Я пытаюсь выполнить линейную регрессию для такой модели:

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

Так что никаких ограничений на б.

Я знаю, что пакет glmc можно использовать для применения ограничений, но я не смог определить, как применить его для моих ограничений. Я также знаю, что contr.sum можно использовать, например, чтобы все коэффициенты суммировались с 0, но я не хочу этого делать. solve.QP() кажется другой возможностью, где установка meq=0 можно использовать так, чтобы все коэффициенты были>=0 (опять же, здесь не моя цель).

Примечание. Решение должно иметь возможность обрабатывать значения NA в векторе ответов Y, например, с помощью:

Y <- runif(100, -1.0, 1.0)
Y[c(2,5,17,56,37,56,34,78)] <- NA

Answer 1

solve.QP могут быть переданы произвольные линейные ограничения, так что, безусловно, может использоваться для моделирования ваших ограничений a+c >= 0 а также c >= 0,

Во-первых, мы можем добавить столбец 1 к X чтобы захватить член перехвата, а затем мы можем повторить стандартную линейную регрессию с 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.

Вы можете обрабатывать пропущенные значения в Y или же X2 так же, как lm Функция делает, удаляя оскорбительные наблюдения. Вы можете сделать что-то вроде следующего в качестве шага предварительной обработки:

has.missing <- rowSums(is.na(cbind(Y, X2))) > 0
Y <- Y[!has.missing]
X2 <- X2[!has.missing,]