For K categories the joint multinomial model uses the
softmax probabilities
\pi_{ik}=\frac{\exp(\eta_{ik})}{\sum_{h=1}^K\exp(\eta_{ih})}.
mgcv::multinom(K-1) selects this family. The response
can be a factor or a matrix of category counts. Count rows may contain
multiple trials.
Fit and predict probabilities
set.seed(22)
X <- matrix(rnorm(400*20), 400, 20)
eta <- cbind(1.2*X[, 1], -X[, 1]+X[, 5], -X[, 5])
prob <- exp(eta-apply(eta, 1, max)); prob <- prob/rowSums(prob)
y <- factor(vapply(seq_len(nrow(X)), function(i)
sample.int(3, 1, prob=prob[i, ]), integer(1)), levels=1:3)
fit <- L0glm.fit(X, y, family=mgcv::multinom(K=2),
lambda="gic", grouped=TRUE, relaxed=TRUE, maxit=150, tol=1e-8)
predict(fit, X[1:5, ], type="response")
#> 1 2 3
#> [1,] 0.1912 0.52463 0.2842
#> [2,] 0.9579 0.00259 0.0395
#> [3,] 0.7564 0.12808 0.1155
#> [4,] 0.4121 0.23701 0.3509
#> [5,] 0.0574 0.92558 0.0170Interpret CLR coefficients
Adding the same constant to every class predictor leaves softmax probabilities unchanged. L0glm reports centered-log-ratio (CLR) coefficients, with each row summing to zero. A positive coefficient increases that category’s log probability relative to the geometric mean of category probabilities.
coef(fit)[1:6, , drop=FALSE]
#> 6 x 3 sparse Matrix of class "dgCMatrix"
#> 1 2 3
#> Intercept 0.00987 -0.216 0.2066
#> 1.24312 -1.175 -0.0677
#> . . .
#> . . .
#> . . .
#> -0.10846 1.178 -1.0696
max(abs(rowSums(as.matrix(coef(fit)))))
#> [1] 1.8e-16grouped=TRUE selects whole predictor rows using
class-symmetric row norms; the model is fitted jointly.
grouped=FALSE promotes sparsity in individual displayed CLR
entries, subject to their zero-sum constraint. Grouped relaxed refitting
admits all class coefficients in selected rows; ungrouped refitting
retains the selected CLR-entry pattern.
Explicit starting matrices are CLR-centered. Rows that are zero in
every channel are excluded in advance, with zero rows restored in the
original output dimensions; mandatory rows remain eligible. Scalar
start=1 centers to zero in CLR space and does not lock out
every predictor.
Tune the grouped model
gic <- L0glm.tune(X, y, family=mgcv::multinom(2), tune="gic",
nlambda=10, grouped=TRUE, relaxed=TRUE, maxit=150)
cv <- L0glm.tune(X, y, family=mgcv::multinom(2), tune="cv",
nfolds=5, nlambda=10, grouped=TRUE, relaxed=TRUE, maxit=150)
predict(cv$fit_1se, X[1:5, ], type="response")The multinomial
reference gives the joint IRLS equations and parameter counts. See
also mgcv::multinom()
for the family constructor.
