# Restricting model parameters in logistic models in R

Is there any function in R that can solve the problem like this example from the SAS website:

Beginning in SAS 9.3, PROC FMM can be used as an alternative to the LOGISTIC and GENMOD procedures for fitting generalized linear models such as logistic and poisson models. You can fit the model in PROC FMM and use its RESTRICT statement to impose equality or inequality constraints on the model parameters.

For example, in the following logistic model suppose you want to constrain the parameters for X1 and X2 to be equal.

 proc logistic;
model y = x1 x2 x3 x4;
run;


The following statements fit the model in PROC FMM and impose the restriction.

 proc fmm;
model y = x1 x2 x3 x4 / dist=binary link=logit;
restrict x1 1 x2 -1;
run;


To restrict the parameter on X1 to exceed that of X2, use the following RESTRICT statement.

 restrict x1 1 x2 -1 > 0;

• Equality constraints can be accomodated with reparameterization. Equality to a constant can be done with an offset in straight GLM. – Glen_b Apr 26 '13 at 4:06
• For non-negativity, this answer points out it's covered in MASS (the book, rather than the R package), and can be achieved by using box constraints with a suitable optimizer passed to glm.fit. With reparameterization that can also cover many cases like $\theta_1 - \theta_2 \geq 0$ – Glen_b Apr 26 '13 at 4:13
• I'm glad. I didn't really think those pointers would be much help. – Glen_b Apr 26 '13 at 5:54

Sorry to keep you waiting for over 5 years, but I just came up with a cool little hack that apparently solves for this. It's based on iteratively updating a constant offset until convergence. The R code below simulates data from a logistic model with $$logit(\Pr(Y= 1 | x_1, x_2)) = -2 + .5 x_1 + .5x_2$$, and creates a situation where both regression parameters are equal.

set.seed(124)
N <- 500
sim_df <- data.frame(x1 = rnorm(N))
sim_df$$x2 <- 1.3 + .5 + sim_df$$x1 + rnorm(N)
cor(sim_df$$x1, sim_df$$x2) # to keep it interesting

sim_df$$z <- with(sim_df, -2 + .5 * x1 + .5 * x2) sim_df$$u <- runif(N)

sim_df$$y <- with(sim_df, as.numeric(u < plogis(z))) table(sim_df$$y)

my_model <- glm(y ~ x1 + x2, data = sim_df, family = "binomial")
summary(my_model)


Even though the true values of the coefficients are the same, the estimated values are clearly different:

Call:
glm(formula = y ~ x1 + x2, family = "binomial", data = sim_df)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-1.8887  -0.8124  -0.5640   0.9315   2.5558

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -1.8947     0.2382  -7.956 1.78e-15 ***
x1            0.3558     0.1557   2.284   0.0224 *
x2            0.5467     0.1109   4.928 8.30e-07 ***
---


Let's initialize based on the above model and iterate:

beta1 <- coef(my_model)[["x1"]]
prev_beta1 <- 0

while (abs(prev_beta1 -  beta1) > .00001) {
my_model <- glm(y ~ x1 + offset(x2 * beta1), data = sim_df,
family = "binomial")
prev_beta1 <- beta1
beta1 <- coef(my_model)[["x1"]]
cat("previous beta1: ", prev_beta1, "new beta1: ", beta1, "\n")
}

  -- This is output:
previous beta1:  0.3557504 new beta1:  0.5203537
previous beta1:  0.5203537 new beta1:  0.3778784
previous beta1:  0.3778784 new beta1:  0.5007586
previous beta1:  0.5007586 new beta1:  0.3944498
previous beta1:  0.3944498 new beta1:  0.4861758
previous beta1:  0.4861758 new beta1:  0.406849
previous beta1:  0.406849 new beta1:  0.4753156

summary(my_model)

Call:
glm(formula = y ~ x1 + offset(x2 * beta1), family = "binomial",
data = sim_df)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-1.8289  -0.8224  -0.5777   0.9677   2.4802

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -1.6990     0.1077  -15.78  < 2e-16 ***
x1            0.4435     0.1239    3.58 0.000343 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 541.75  on 499  degrees of freedom
Residual deviance: 527.63  on 498  degrees of freedom
AIC: 531.63

Number of Fisher Scoring iterations: 4


Note that the final model with the offset doesn't include $$x_2$$, but you can see below that it's in the model, with the same coefficient value as $$x_1$$:

# compare:
predict(my_model)[1]
# vs
coef(my_model)[1] +
coef(my_model)[["x1"]] * sim_df[1, "x1"] +
coef(my_model)[["x1"]] * sim_df[1, "x2"]