Restricting model parameters in logistic models in R

Question

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;

Answer 1

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"]