Как ограничить коэффициенты регрессии Пуассона положительными в R?

Question

Предположим, что следующая ситуация:
У меня есть данные переменной Y, которую я предполагаю распределить по Пуассону. У меня также есть данные переменной X за тот же период времени, и каждое наблюдение представляет собой определенное событие. Я полагаю, что значения Y исходят из двух разных воздействий, поэтому я разделяю каждое наблюдение Y_i на два распределенных по Пуассону Y_i1 и Y_i2, но у меня все еще есть только наблюдения относительно общего количества Y_i. Я также предполагаю, что события (представленные X) оказывают долгосрочное влияние на Y_i1, и у меня есть оценки параметра lambda_i2.

Таким образом, моя формула регрессии -
fml=Y_i ~ b_1*X_i+....+b_n*X_(i-n+1) + offset(lambda_i2) -1 с n> = 24.
Это означает, что последние 24 (или более, из-за долгосрочного эффекта) значения X влияют на значение Y_i1 аддитивным способом, и у меня нет перехвата (b_0 = 0).

Я сделал матрицу m, строки которой представляют Y_i, все ее 24 (или более) регрессоры для каждого наблюдения Y_i и соответствующей оценки lambda_i2.

Теперь я использовал glm(fml, family=poisson(link="identity"), data=m) и попробовал его для разных значений n (= 24,48,36, ...).

Всегда, некоторые из коэффициентов получили отрицательные значения, которые не имеют смысла в интерпретации. (События, представленные X может иметь только положительный или NONE эффект на величину Y.)

Это приводит к моему вопросу:

Как я могу использовать ограничение b_i >=0 в моей модели?

В моем предыдущем исследовании я нашел функцию glmc(), но я не уверен, как включить мое ограничение здесь.

В качестве альтернативы я также подумал об анализе этой модели байесовским способом, но пока не нашел байесовскую версию glm() для распределения Пуассона, так что я могу указать предыдущий для b_i самостоятельно , (Тогда я мог бы включить положительность в предыдущем.)

Есть ли у вас идеи?

Это выдержка из моих данных и мой кода:

y=c(279,623, 1025, 1701, 1862, 2544, 2308, 2231, 2234, 2550, 2698, 2805, 3510, 3032, 2746, 2074, 1062, 513, 226, 116, 87, 79, 116, 335, 594, 1081, 1425, 1775, 2056, 2387, 2337, 2354, 2665, 2406, 2433, 2550, 2820, 3655, 4566, 2330, 1267, 531, 280, 148, 92, 89, 141, 458, 852, 1214) 
X=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.88349136, 0.54951680, 0.13306912, 0.15321180, 0.00000000, 1.42569128, 0.55808054, 0.65486418, 0.27530564, 0.24813572, 0, 0, 2.09889028, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.18947898, 0.17347032, 0.94538886, 0.03334654, 0.05593732, 0.00000000, 0.99772264, 0.11121918, 0, 1.41673120, 0.27375384, 0, 0, 0, 0, 5.67487576, 0, 0, 0, 0, 0, 0, 0, 0, 1.55642510, 0.98419866, 0.50992652) 
lambda=c(253.5, 562.5, 1053.0, 1645.0, 2064.5, 2215.0, 2503.0, 2443.0, 2514.5, 2701.0, 2972.5, 3035.5, 3422.5, 3295.0, 2882.5, 2094.0, 1211.0, 579.5, 265.5, 155.0, 112.5, 82.5, 117.5, 306.0, 627.0, 1021.0, 1463.5, 1722.5, 2017.5, 2146.5, 2209.0, 2231.5, 2265.0, 2320.0, 2442.0, 2507.0, 2957.0, 3674.0, 3345.5, 2285.0, 1265.0, 555.5, 252.0, 145.5, 86.5, 90.5, 148.0, 362.0, 738.0, 1137.5) 

regressors=function(n,x){ 
    m=length(x)-n+1; 
    r=matrix(0,m,n); 
    for (i in 0:(n-1)){ r[,(i+1)]=x[(n-i):(length(x)-i)]} 
    return(r); 
} 
r=regressors(24,X); 
reg=cbind(y,data.frame(r),lambda); 
fml=as.formula(paste("y~", paste(colnames(reg)[2:25], collapse = "+"), "+offset(lambda)-1")); 
g=glm(fml, poisson(link="identity"), data=reg); %this leads to negative coefficients 
obj=function(b){-sum(dpois(y, r%*%b, log=TRUE))} 
st=coef(lm(fml, data=reg)); 
opt=optim(st, obj); % this is where the error occurs 

регрессоров() является функцией я написал для вычисления регрессоров (это приводит к матрице г с п столбцами и 50 строк и каждой строкой i представляет регрессоры y_i).

Answer 1

Попробуйте первые принципы.

# generate random input data 
set.seed(123) 
n <- 100 
x <- 1:n 
X <- cbind(1, x) 
b <- c(0.1, 3) 
y <- rpois(n, X %*% b) 

# log likelihood objective function 
obj <- function(b) -sum(dpois(y, X %*% b, log = TRUE)) 

# as a check try with no constraints - these two should give the same coefs 
glm(y ~ X + 0, family = poisson("identity")) 
st <- coef(lm(y ~ X + 0)); optim(st, obj) 

# now add lower bounds ensuring starting value is in feasible region 
optim(pmax(st, 1), obj, lower = c(0, 0), method = "L-BFGS-B") 

Примечание 1: Будьте осторожны, если ваши оценки параметров на границе допустимой области, как это происходит в стесненном примере выше.

Примечание 2: Ниже приведена переделка кода, добавленного к вопросу. y, X, lambda и regressors как в вопросе. Замечание: добавление ограничения приводит к тому, что многие коэффициенты равны нулю.

r <- regressors(24,X) 
reg <- cbind(y,data.frame(r),lambda) 

fml <- y ~ . - lambda + offset(lambda) 

# check that g and opt give the same coefs 

g <- glm(fml, poisson(link = "identity"), data = reg) 

obj <- function(b)-sum(dpois(y, r%*%b + lambda, log = TRUE)) 
st <- coef(lm(fml, data=reg)) 
opt <- optim(st, obj, method = "BFGS", control = list(maxit = 500)) 

optim(pmax(coef(g), 1), obj, method = "L-BFGS-B", lower = 0 * st)