fitting a regression with coefficients constrained to sum to 1

softmax_regression.R

## https://fediscience.org/@[email protected]/113439584737810431

library(MASS)   # for eqscplot
library(RTMB)

set.seed(101)
n <- 1000
X <- cbind(
  x1 = rnorm(n , 10, 2),
  x2 = rnorm(n, 7, 1),
  x3 = rnorm(n, 15, 3),
  x4 = runif(n, 5, 20),
  x5 = rnorm(n, 8, 1.5)
)
# note true data generating process coefficients add up
# to more than 1 and have a negative, to make it interesting
true_coefs <- c(0.2, 0.6, -0.1, 0, 0.4)

y <- X %*% true_coefs + rnorm(n, 0, 2)

softmax <- function(x) {
    v <- c(0, x)
    ee <- exp(v)
    ee/sum(ee)
}

nllfun <- function(pars) {
    getAll(pars)
    beta <- softmax(beta0)
    y %~% dnorm(X %*% beta, exp(logsd))
}
pars0 <- list(beta0 = rep(0, ncol(X)), logsd = 0)
ff <- MakeADFun(nllfun, pars0, silent = TRUE)

fit <- with(ff, nlminb(unlist(pars0), fn, gr))
beta_fit <- softmax(fit$par[1:ncol(X)])

f <- function(newx, main = ""){
  eqscplot(true_coefs, newx, bty = "l", 
           xlab = "Coefficients from original data generating process",
           ylab = "Coefficients from model",
           main = main) 
  abline(0, 1, col = "grey")
}

f(beta_fit)

oops, there may be a small indexing issue here. I think lines 32, 36 should be ncol(X)-1? Not sure why this didn't trigger errors ...