Automatic Response Category Combination in Multinomial Logistic Regression

TL;DR

This paper introduces a penalized likelihood approach for multinomial logistic regression that automatically combines response categories, improving model interpretability and prediction accuracy.

Contribution

It develops a novel non-differentiable penalty and an ADMM algorithm for efficient estimation and category combination in multinomial logistic regression.

Findings

Effective category combination achieved in simulations

Algorithm converges reliably in practice

Improved prediction performance over traditional methods

Abstract

We propose a penalized likelihood method that simultaneously fits the multinomial logistic regression model and combines subsets of the response categories. The penalty is non differentiable when pairs of columns in the optimization variable are equal. This encourages pairwise equality of these columns in the estimator, which corresponds to response category combination. We use an alternating direction method of multipliers algorithm to compute the estimator and we discuss the algorithm's convergence. Prediction and model selection are also addressed.

Tables

Table 1. Table 1: Results of simulation using setting 1 when ℒ = 𝒮 ℒ 𝒮 \operatorname{\mathcal{L}}=\operatorname{\mathcal{S}} . Standard errors are presented in parenthesis. The columns Oracle, VL, EN, and GMR report to the average KL divergence for oracle tuned GFMR, GFMR using validation likelihood, EN, and GMR methods respectively. The column VL-Oracle reports the average difference for GFMR using validation likelihood and oracle tuned GFMR. Similar notation is used for comparisons between VL and EN, and VL and GMR. The columns labeled Oracle Groups, and VL Groups, contain the average number of unique regression coefficient vectors when using oracle tuned GFMR, and GFMR using validation likelihood.

	$\delta$	Oracle	VL	VL-Oracle	Oracle Groups	VL Groups	EN	VL-EN	GMR	VL-GMR
$n=50$	0.10	0.011 (0.001)	0.023 (0.000)	0.011 (0.002)	1.95	2.98	0.050 (0.003)	-0.027 (0.003)	0.049 (0.003)	-0.027 (0.003)
	0.25	0.057 (0.002)	0.076 (0.002)	0.019 (0.002)	2.98	3.91	0.094 (0.003)	-0.017 (0.004)	0.091 (0.003)	-0.015 (0.003)
	0.50	0.137 (0.003)	0.164 (0.005)	0.027 (0.003)	3.51	4	0.191 (0.005)	-0.027 (0.004)	0.188 (0.004)	-0.024 (0.004)
	1.00	0.259 (0.005)	0.292 (0.005)	0.033 (0.004)	4	4	0.399 (0.010)	-0.107 (0.009)	0.399 (0.010)	-0.106 (0.008)
$n=100$	0.10	0.011 (0.001)	0.019 (0.001)	0.008 (0.001)	2.27	3.36	0.033 (0.002)	-0.013 (0.001)	0.032 (0.002)	-0.012 (0.001)
	0.25	0.050 (0.001)	0.061 (0.001)	0.008 (0.001)	2.27	3.36	0.033 (0.002)	-0.009 (0.001)	0.069 (0.002)	-0.012 (0.001)
	0.50	0.098 (0.002)	0.114 (0.004)	0.016 (0.003)	3.86	4	0.141 (0.005)	-0.027 (0.003)	0.142 (0.004)	-0.027 (0.003)
	1.00	0.150 (0.004)	0.164 (0.004)	0.015 (0.002)	4	4	0.222 (0.006)	-0.058 (0.005)	0.225 (0.006)	-0.061 (0.005)

Table 2. Table 2: Results of simulation using setting 1 when ℒ = { ( 1 , 2 ) , ( 2 , 3 ) , ( 3 , 4 ) } ℒ 1 2 2 3 3 4 \operatorname{\mathcal{L}}=\{(1,2),(2,3),(3,4)\} . The columns Oracle, and VL report to the average KL divergence for oracle tuned GFMR and GFMR using validation likelihood methods respectively. The column VL-Oracle reports the average difference for GFMR using validation likelihood and oracle tuned GFMR. Similar notation is used for comparison VL and GMR. The columns labeled Oracle Groups, and VL Groups, contain the average number of unique regression coefficient vectors when using oracle tuned GFMR, and GFMR using validation likelihood.

	$\delta$	Oracle	VL	VL-Oracle	Oracle Groups	VL Groups	VL-GMR
$n=50$	0.10	0.0121 (0.000)	.025 (0.003)	0.013 (0.002)	3.1	3.1	-0.024 (0.003)
	0.25	0.065 (0.001)	0.0.086 (0.003)	0.028 (0.003)	3.4	3.25	-0.006 (0.003)
	0.50	0.156 (0.004)	0.198 (0.005)	0.040 (0.005)	3.95	3.62	0.009 (0.004)
	1.00	0.189 (0.003)	0.245 (0.003)	0.054 (0.001)	4	3.77	-0.017 (0.010)
$n=100$	0.10	0.012 (0.000)	0.020 (0.002)	0.008 (0.001)	2.91	3.15	-0.011 (0.001)
	0.25	0.058 (0.001)	0.069 (0.001)	0.011 (0.002)	3.65	3.41	0.000 (0.002)
	0.50	0.107 (0.003)	0.127 (0.005)	0.020 (0.004)	3.99	3.86	-0.014 (0.003)
	1.00	0.156 (0.004)	0.191 (0.005)	0.035 (0.004)	4	3.99	-0.034 (0.003)

Table 3. Table 3: Number of replications out of 100 that produced the number of unique regression coefficient vectors in setting 1. Ordered penalty set references that the penalty set used is 𝒮 𝒮 \operatorname{\mathcal{S}} , while unordered penalty set references ℒ = { ( 1 , 2 ) , ( 2 , 3 ) , ( 3 , 4 ) } ℒ 1 2 2 3 3 4 \operatorname{\mathcal{L}}=\{(1,2),(2,3),(3,4)\} .

	$\delta$	Penalty Set	1 Vector	2 Vectors	3 Vectors	4 Vectors
$n=50$	0.10	Ordered	0	6	78	16
		Unordered	64	6	1	29
	0.25	Ordered	0	5	65	30
		Unordered	40	3	3	54
	0.50	Ordered	0	1	36	63
		Unordered	16	0	1	83
	1.00	Ordered	0	0	12	88
		Unordered	0	0	0	100
$n=10$	0.10	Ordered	0	8	68	24
		Unordered	49	12	2	37
	0.25	Ordered	0	1	57	42
		Unordered	24	3	4	69
	0.50	Ordered	0	0	14	86
		Unordered	4	1	0	95
	1.00	Ordered	0	0	1	99
		Unordered	0	0	0	100

Table 4. Table 4: Results of simulation using setting 2 when ℒ = 𝒮 ℒ 𝒮 \operatorname{\mathcal{L}}=\operatorname{\mathcal{S}} . Standard errors are presented in parenthesis. The columns Oracle, VL, EN, and GMR report to the average KL divergence for oracle tuned GFMR, GFMR using validation likelihood, EN, and GMR methods respectively. The column VL-Oracle reports the average difference for GFMR using validation likelihood and oracle tuned GFMR. Similar notation is used for comparisons between VL and EN, and VL and GMR. The columns labeled Oracle Groups, and VL Groups, contain the average number of unique regression coefficient vectors when using oracle tuned GFMR, and GFMR using validation likelihood.

	$\delta$	Oracle	VL	VL-Oracle	Oracle Groups	VL Groups	EN	VL-EN	GMR	VL-GMR
$n=50$	0.10	0.008 (0.000)	0.023 (0.003)	0.014 (0.003)	1.99	2.81	0.044 (0.003)	-0.022 (0.002)	0.043 (0.003)	-0.023 (0.003)
	0.25	0.043 (0.001)	0.056 (0.002)	0.013 (0.002)	2.45	3.75	0.083 (0.003)	-0.026 (0.003)	0.081 (0.003)	-0.025 (0.003)
	0.50	0.118 (0.003)	0.146 (0.003)	0.026 (0.002)	3.07	4	0.169 (0.003)	-0.023 (0.004)	0.167 (0.003)	-0.021 (0.003)
	1.00	0.225 (0.006)	0.243 (0.007)	0.017 (0.004)	4	4	0.364 (0.008)	-0.120 (0.008)	0.372 (0.009)	-0.120 (0.008)
$n=100$	0.10	0.008 (0.000)	0.014 (0.001)	0.006 (0.001)	2.04	3.39	0.027 (0.001)	-0.012 (0.001)	0.027 (0.002)	-0.012 (0.001)
	0.25	0.037 (0.001)	0.048 (0.001)	0.010 (0.001)	2.90	3.91	0.060 (0.002)	-0.012 (0.002)	0.059 (0.002)	-0.012 (0.002)
	0.50	0.076 (0.002)	0.089 (0.003)	0.012 (0.002)	3.85	4	0.116 (0.002)	-0.027 (0.002)	0.115 (0.003)	-0.026 (0.002)
	1.00	0.124 (0.004)	0.131 (0.004)	0.006 (0.001)	4	4	0.201 (0.007)	-0.070 (0.005)	0.204 (0.006)	-0.073 (0.005)

Table 5. Table 5: Results of simulation using setting 2 when ℒ = { ( 1 , 2 ) , ( 2 , 3 ) , ( 3 , 4 ) } ℒ 1 2 2 3 3 4 \operatorname{\mathcal{L}}=\{(1,2),(2,3),(3,4)\} . The columns Oracle, and VL report to the average KL divergence for oracle tuned GFMR and GFMR using validation likelihood methods respectively. The column VL-Oracle reports the average difference for GFMR using validation likelihood and oracle tuned GFMR. Similar notation is used for comparison VL and GMR. The columns labeled Oracle Groups, and VL Groups, contain the average number of unique regression coefficient vectors when using oracle tuned GFMR, and GFMR using validation likelihood.

	$\delta$	Oracle	VL	VL-Oracle	Oracle Groups	VL Groups	VL-GMR
$n=50$	0.10	0.008 (0.003)	0.023 (0.004)	0.015 (0.003)	2.87	2.11	-0.025 (0.003)
	0.25	0.040 (0.000)	0.057 (0.002)	0.016 (0.001)	3.17	3.56	-0.024 (0.003)
	0.50	0.106 (0.003)	0.135 (0.004)	0.029 (0.003)	3.79	3.56	-0.032 (0.004)
	1.00	0.189 (0.003)	0.245 (0.003)	0.054 (0.001)	4	3.77	-0.128 (0.010)
$n=100$	0.10	0.008 (0.000)	0.015 (0.001)	0.007 (0.001)	2.91	3.15	-0.011 (0.001)
	0.25	0.035 (0.001)	0.046 (0.002)	0.011 (0.002)	3.43	3.35	-0.013 (0.002)
	0.50	0.064 (0.003)	0.073 (0.003)	0.009 (0.001)	3.89	3.75	-0.042 (0.002)
	1	0.101 (0.003)	0.119 (0.004)	0.018 (0.002)	4	3.89	-0.085 (0.006)

Table 6. Table 6: Number of replications out of 100 that produced the number of unique regression coefficient vectors in setting 1. Ordered penalty set references that the penalty set used is 𝒮 𝒮 \operatorname{\mathcal{S}} , while unordered penalty set references ℒ = { ( 1 , 2 ) , ( 2 , 3 ) , ( 3 , 4 ) } ℒ 1 2 2 3 3 4 \operatorname{\mathcal{L}}=\{(1,2),(2,3),(3,4)\} .

	$\delta$	Penalty Set	1 Vector	2 Vectors	3 Vectors	4 Vectors
$n=50$	0.10	Ordered	0	2	83	14
		Unordered	63	4	4	29
	0.25	Ordered	0	4	69	27
		Unordered	49	2	4	45
	0.50	Ordered	0	2	40	58
		Unordered	29	3	0	68
	1.00	Ordered	0	0	22	78
		Unordered	0	0	0	100
$n=100$	0.10	Ordered	0	7	71	22
		Unordered	55	14	3	28
	0.25	Ordered	0	3	59	38
		Unordered	33	4	3	60
	0.50	Ordered	0	0	25	75
		Unordered	4	1	1	94
	1.00	Ordered	0	0	11	89
		Unordered	0	0	0	100

Table 7. Table 7: The fraction of the 100 replications specific group structures are selected for each N 𝑁 N , δ 𝛿 \delta combination using the two-step method . The label One-Step indicates that the correct group structure is still a possibility if the correct fusion was done with an additional combination.

	$N=50$		$N=75$
	$\delta =1$	$\delta =3$	$\delta =1$	$\delta =3$
1 Group	19/100	0/100	3/100	0/100
2 Groups (Correct)	58/100	71/100	80/100	83/100
2 Groups (Incorrect)	0/100	0/100	0/100	0/100
3 Groups (One-Step)	20/100	21/100	16/100	15/100
3 Groups (Incorrect)	0/100	0/100	0/100	0/100
4 Groups	3/100	8/100	1/100	2/100

Table 8. Table 8: The fraction of the 100 replications specific group structures are selected for each N 𝑁 N , δ 𝛿 \delta combination when all possible response category combinations were used as candidate models.

	$N=50$		$N=75$
	$\delta =1$	$\delta =3$	$\delta =1$	$\delta =3$
1 Group	0/100	0/100	0/100	0/100
2 Groups (Correct)	46/100	59/100	62/100	82/100
2 Groups (Incorrect)	25/100	4/100	14/100	0/100
3 Groups (One-Step)	28/100	29/100	24/100	17/100
3 Groups (Incorrect)	0/100	0/100	0/100	0/100
4 Groups	1/100	8/100	0/100	1/100

Table 9. Table 9: The response categories found by GFMR with tuning parameter selection using validation likelihood on the 1996 United States election data.

Group	$ℒ=𝒮$	$ℒ=\{(1,2),\dots ,(6,7)\}$
1	Strong Republican	Strong Republican
2	Weak Republican	Weak Republican
3	Independent Republican	Independent Republican Independent Independent Democrat
4	Independent	Weak Democrat
5	Independent Democrat	Strong Democrat
6	Weak Democrat
7	Strong Democrat

Table 10. Table 10: The response categories found by GFMR using the two step method for model selection on the 1996 United States election data.

Group	$ℒ=𝒮$	$ℒ=\{(1,2),\dots ,(6,7)\}$
1	Strong Republican Weak Republican Independent Republican Independent Democrat	Strong Republican Weak Republican
2	Independent	Independent Republican Independent Independent Democrat
3	Strong Democrat Weak Democrat	Strong Democrat Weak Democrat

Equations

$$Y_i\sim \mathrm{Multinomial}\left(n_i,{\pi }_{i1}^{\ast },\dots ,{\pi }_{iC}^{\ast }\right),i=1,\dots ,n$$

$${\pi }_{ij}^{\ast }=\frac{\exp (x_i^{\prime }{\beta }_j^{\ast })}{{\sum }_{m=1}^C\exp (x_i^{\prime }{\beta }_m^{\ast })},j\in C,$$

$$\ell ({\beta }_1,\dots ,{\beta }_{C-1})=\sum _{i=1}^n\sum _{j=1}^{C-1}{y}_{ij}{x}_i^{\prime }{\beta }_j-\sum _{i=1}^n{n}_i\log \left\{1+\sum _{r=1}^{C-1}\exp \left(x_i^{\prime }{\beta }_r\right)\right\}.$$

$$({\hat{\beta }}_1,\dots ,{\hat{\beta }}_C)=\mathrm{argmin}{\text{displaylimits}}_{({\beta }_1,\dots ,{\beta }_C)\in {\mathbb{R}}^p\times \cdots \times {\mathbb{R}}^p}$$

$$\mathrm{subject}\mathrm{to}{\beta }_C=\vec{0},$$

$$Q(\lambda )=\sum _{k=1}^K\sum _{i\in V_k}\sum _{l=1}^C{y}_{il}\log \left({\hat{\pi }}_l^{(-V_k)}(x_i,\lambda )\right),$$

$$\mathrm{argmin}{\text{displaylimits}}_{{\beta }_j\in {\mathbb{R}}^p,j=1,\dots ,C,Z_{j,m}\in {\mathbb{R}}^p,(j,m)\in L}$$

$$\text{mbox}subjectto{\beta }_C=\vec{0}\text{mbox}and{Z}_{j,m}={\beta }_j-{\beta }_m,(j,m)\in L.$$

$$L_{\rho }(\beta ,Z,W)=$$

$$+\sum _{(j,m)\in L}\left\{\lambda |Z_{j,m}{|}_2+W_{j,m}^T({\beta }_j-{\beta }_m-Z_{j,m})+\frac{\rho }{2}|{\beta }_j-{\beta }_m-Z_{j,m}{|}_2^2\right\},$$

$${\beta }^{(k+1)}$$

$$Z^{(k+1)}$$

$$W_{j,m}^{(k+1)}$$

$$Z_{j,m}^{(k+1)}=\left({\beta }_j^{(k+1)}-{\beta }_m^{(k+1)}+\frac{1}{\rho }{W}_{j,m}^{(k)}\right){\left(1-\frac{\lambda }{\rho |{\beta }_j^{(k+1)}-{\beta }_m^{(k+1)}+\frac{1}{\rho }{W}_{j,m}^{(k)}{|}_2}\right)}_{+},(j,m)\in L,$$

$${\hat{\beta }}_j=\frac{{\sum }_{m\in C\setminus \{j\}}{\beta }_m^{(K)}I(Z_{j,m}^{(K)}=\vec{0})+{\beta }_j^{(K)}}{{\sum }_{m\in C\setminus \{j\}}I(Z_{j,m}^{(K)}=\vec{0})+1},j=1,\dots ,C-1,$$

$$KL(\hat{\pi },{\pi }^{\ast })=\frac{1}{1000}\sum _{i=1}^{1000}\sum _{k=1}^4\log \left(\frac{{\hat{\pi }}_k(x_i,\lambda )}{{\pi }_k^{\ast }(x_i)}\right){\hat{\pi }}_k(x_i,\lambda ),$$

Full text

Key Words: Multinomial logistic regression; Fusion penalty; Response category reduction.

1 Introduction

We propose a new way to fit the multinomial logistic regression model. Let $x_{i}=(1,x_{i2},\ldots,x_{ip})^{\prime}\in\mathbb{R}^{p}$ be the non-random values of the predictors for the $i$th subject and let $y_{i}=(y_{i1},\ldots,y_{iC})^{\prime}\in\mathbb{R}^{C}$ be the observed response category counts for the $i$th subject $(i=1,\ldots,n)$. The model assumes that $y_{i}$ is a realization of the random vector

[TABLE]

where $n_{i}$ is the index of the multinomial experiment for the $i$th subject and $\pi^{*}_{ij}$ is the unknown probability that the response is category $j$ for the $i$th subject $(i,j)\in\{1,\ldots,n\}\times\operatorname{\mathcal{C}}$, where $\operatorname{\mathcal{C}}=\{1,2,\ldots,C\}$. The model also assumes that $Y_{1},\ldots,Y_{n}$ are independent. Using the baseline category parameterization,

[TABLE]

where $\beta_{1}^{*},\ldots,\beta_{C-1}^{*}$ are unknown regression coefficient vectors and $\beta_{C}^{*}=\vec{0}$. Other constraints could be used to make $\beta^{*}_{1},\ldots,\beta_{C}^{*}$ identifiable, e.g. Zhu and Hastie (2004) used $\sum_{j=1}^{C}\beta^{*}_{j}=0$.

Let $\ell:\mathbb{R}^{p}\times\cdots\times\mathbb{R}^{p}\rightarrow\mathbb{R}$ be the log likelihood with additive terms that do not contain parameters dropped:

[TABLE]

There are a total of $p(C-1)$ unknown parameters. For more information on the model see Agresti (2012, Chapter 8).

Shrinkage or regularized estimation is natural in this setting when $p(C-1)$ is large. Zhu and Hastie (2004) proposed ridge-penalized likelihood estimation, and Vincent and Hansen (2014) proposed sparse group-lasso penalized likelihood estimation that encourages estimates of the matrix $(\beta_{1},\ldots,\beta_{C-1})\in\mathbb{R}^{p\times C-1}$ that have 0’s in all entries of some rows. These rows correspond to explanatory variables that are estimated to be irrelevant. A similar procedure was studied by Simon et al. (2013).

Interpreting the regression coefficients when the response has three or more categories is difficult because the coefficient values depend on the baseline category choice (or depend on another constraint employed to make the parameters identifiable). Interpretation in the special case of binomial logistic regression ($C=2$) is much simpler. In some applications, we can create a two-category response by subjectively grouping the original response categories, e.g. if the $j$th response category was of interest, one could create two groups of categories $\{j\}$ and $\operatorname{\mathcal{C}}\setminus\{j\}$, and perform binomial logistic regression. Without subjective information, we can view response category grouping as a model selection problem, where the selected submodel with grouped categories fits the counts in our data nearly as well as the full model.

We propose a penalized likelihood procedure that uses fusion penalties to encourage fitted models with grouped response categories. Model selection is addressed using both $K$-fold cross validation and AIC. Through simulation, we show that in settings where the true response categories are grouped, our method using cross validation for tuning parameter selection performs better than competitors at predicting the true response category probabilities. We also show that in certain settings, our method using AIC for model selection excels at selecting the true grouped response categories.

2 Method

We propose the penalized likelihood estimates defined by

[TABLE]

where $\mathcal{L}$ is a user-selected subset of $\operatorname{\mathcal{S}}=\{(a,b)\in\operatorname{\mathcal{C}}\times\operatorname{\mathcal{C}}:a>b\}$, $\lambda\geq 0$ is a tuning parameter, and $|\cdot|_{2}$ is the vector 2-norm. Since the penalty only depends on differences of $\beta$ vectors, the estimated response category probabilities are invariant to the choice of the baseline category for a given $\operatorname{\mathcal{L}}$.

To exploit a natural ordering (if any) of the response categories, one could set $\operatorname{\mathcal{L}}=\{(1,2),(2,3),\ldots,(C-1,C)\}$. Without prior information about the similarity between response categories, one could set $\operatorname{\mathcal{L}}=\operatorname{\mathcal{S}}$.

Suppose that $\operatorname{\mathcal{L}}=\operatorname{\mathcal{S}}$. The objective function in (4) is non-differentiable $(\beta_{1},\ldots,\beta_{C})$ has at least one pair of equal components, $\beta_{m}=\beta_{j}$, $m\neq j$. This encourages an increasing number of pairs of equal element vectors in $(\hat{\beta}_{1},\ldots,\hat{\beta}_{C})$ as $\lambda$ increases. If $\hat{\beta}_{m}=\hat{\beta}_{j}$, then the corresponding estimated response category probabilities are also equal. In effect, response categories $j$ and $m$ are combined.

If $\lambda$ is sufficiently large, then $\hat{\beta}_{1}=\cdots=\hat{\beta}_{C}=0$. In this uninteresting edge-case, all $\hat{\pi}_{ij}$’s are equal to $1/C$.

We select $\lambda$ using $K$-fold cross validation maximizing the validation log-likelihood. Let $\operatorname{\mathcal{V}}_{k}$ be the set of indices in the $k$-th fold. Specifically we maximize the function $Q$, defined by

[TABLE]

where $y_{il}$ is the observed count for response category $l$ for validation observation $i$, and where $\widehat{\pi}^{(-\operatorname{\mathcal{V}}_{q})}_{l}(x_{i},\lambda)$ is the estimated probability validation observation $i$, with predictor values $x_{i}$, has response category $l$ using estimates produced from the training set that omits data with indices in $\operatorname{\mathcal{V}}_{q}$ using tuning parameter $\lambda$. We propose an alternative method of tuning parameter selection using AIC in Section 5.

Our penalty builds on the ideas behind the group lasso penalty (Bakin, 1999; Yuan and Lin, 2006) and the fused lasso penalty (Tibshirani et al., 2005). Alaiz et al. (2013) developed the grouped fused lasso penalty, but their work only fused adjacent groups and was designed for single-response penalized least-squares regression. Both Hocking et al. (2011) and Chen et al. (2015) investigate similar penalties in the context of clustering.

3 Algorithm

We propose to solve (4) using an alternating direction method of multipliers algorithm (ADMM). Boyd et al. (2010) provide an introduction and description of the ADMM algorithm. The proposed ADMM algorithm solves the optimization in (4) by solving the equivalent constrained optimization:

[TABLE]

The optimization in (6) can be written as a special case of equation 3.1 of Boyd et al. (2010). Define $\beta=(\beta_{1},\ldots,\beta_{C-1})$, $Z=\{Z_{j,m}\}_{(j,m)\in\operatorname{\mathcal{L}}}$ and $W=\{W_{j,m}\}_{(j,m)\in\operatorname{\mathcal{L}}}$, where $W_{j,m}$ is the Lagrange multiplier associated with the constraint on $Z_{j,m}$. The augmented Lagrangian is

[TABLE]

where $\rho>0$ is the augmented Lagrangian parameter. We now apply equations 3.2, 3.3, and 3.4 of Boyd et al. (2010) to obtain the ADMM update equations for solving (6):

[TABLE]

where a superscript of $(k)$ denotes the $k$th iterate.

We propose to solve (7) with blockwise coordinate descent, where $\beta_{1},\ldots,\beta_{C-1}$ are the blocks. In each block update, we use Newton’s method, which converges because these subproblems have strongly convex objective functions. Solving (7) is the most computationally expensive step.

The optimization in (8) decouples into $\operatorname{card}(\operatorname{\mathcal{L}})$ optimization problems that can be solved in parallel with the following closed form solutions:

[TABLE]

where $(a)_{+}=\max(a,0)$. This solution is a special case of solving the group lasso penalized least squares problem (Yuan and Lin, 2006).

One could show convergence of the proposed ADMM algorithm by using the theory developed by Mota et al. (2011). This proof relies on the assumption that the undirected graph defined by the vertex set $\operatorname{\mathcal{C}}$ and the edge set $\operatorname{\mathcal{L}}$, is connected. In practice the convergence tolerance and step size $\rho$ can be determined by equation 3.12 and 3.13 respectively in Boyd et al. (2010).

The final iterate $\beta^{(K)}$ will not have pairs of equal component vectors, but subsets of its pairs of component vectors will be very similar. We define our final estimates by

[TABLE]

and $\hat{\beta}_{C}=\vec{0}$. We set $\hat{\beta}_{1},\ldots,\hat{\beta}_{C-1}$ equal to [math] if the 2-norm of the vector is less than $10^{-8}$. A similar approach was taken by Danaher et al. (2013).

4 Prediction Performance Simulations

4.1 Competitors

We present simulation studies that compare the prediction performance of our proposed method, which we call group fused multinomial logistic regression (GFMR), elastic net penalized multinomial logistic regression (EN), and group penalized multinomial logistic regression (GMR) (Friedman et al., 2008; Simon et al., 2013). Both methods used for comparison are implemented in the glmnet package in R. An R package implementing the GFMR methods is currently under development.

4.2 Data Generating Models

In Sections 4.3 and 4.4, we present simulation results where the data was generated such that $x_{i}=(1,\tilde{x}_{i})^{\prime}$, where $\tilde{x}_{1},\ldots,\tilde{x}_{n}$ are drawn independently from $N_{9}(0,I)$. We set $n_{i}=1$ for $i=1,\ldots,n$ and $C=4$. We consider two settings for $\beta^{*}$. In setting 1, $\beta^{*}_{1}=\beta^{*}_{4}=\vec{0}$ and $\beta_{2}^{*}=\beta_{3}^{*}=\vec{\delta}$. In setting 2, $\beta^{*}_{1}=\beta^{*}_{2}=\beta^{*}_{3}=\vec{\delta}$, and $\beta^{*}_{4}=\vec{0}$. We consider $(n,\delta)\in\{50,100\}\times\{0.1,0.25,0.5,1\}$. The observed responses $y_{i},\dots,y_{n}$ are a realization of $Y_{1},\ldots,Y_{n}$ defined by (1).

Results are based on 100 replications. We consider two choices for $\operatorname{\mathcal{L}}$: $\operatorname{\mathcal{L}}=\operatorname{\mathcal{S}}$ and $\operatorname{\mathcal{L}}=\{(1,2),(2,3),(3,4)\}$. The second choice attempts to exploit a natural ordering of the response categories. In each replication, we measure the prediction performance using the Kullback-Leibler (KL) divergence between the estimated response category probabilities and the true response category probabilities, based on 1000 test observations generated from the same distribution as the training data. The KL Divergence is defined by

[TABLE]

where $\pi^{*}_{k}(x_{i})$, and $\widehat{\pi}_{k}(x_{i},\lambda)$ are the true and estimated response category probabilities that testing observation $i$ has response category $k$ using tuning parameter $\lambda$. For both GFMR methods we report the average number of unique regression coefficient vectors observed in the estimator for the 100 replications. The number of unique regression coefficient vectors estimated by the GFMR methods is defined as the unique vectors in $(\hat{\beta}_{1},\ldots,\hat{\beta}_{C})$.

Tuning parameters for GFMR are selected from a subset of $\{10^{-10},10^{-9.98},\ldots,10^{9.98},10^{10}\}$. We also computed the best overall KL divergence we could obtain using the candidate values for the tuning parameter. We call this the oracle tuned value. The tuning parameters for the competing methods were selected using 5-fold cross validation minimizing the validation deviance. For the both EN and GMR, the first tuning parameter was selected by default methods in the glmnet package in R and the second tuning parameter was selected from the set $\{0,0.01,\ldots,0.99,1\}$.

4.3 Results for Setting 1

In Table 1, we present the average KL divergence and the average number of unique regression coefficient vectors estimated by the GFMR method for the simulation using setting 1 when $\operatorname{\mathcal{L}}=\operatorname{\mathcal{S}}$. On average GFMR using validation likelihood had a lower KL divergence than its competitors. Since there should be 2 unique regression coefficient vectors in the estimator, the results show GFMR using oracle tuning on average is overselecting the number of unique coefficient vectors when $\delta=0.25,0.5$, and $1$ for both values of $n$. A similar result occurred for GFMR using validation likelihood for all values of $\delta$ and $n$. This may indicate that we need a different method to select the tuning parameter for GFMR if our interest is in model selection, when $\operatorname{\mathcal{L}}=\operatorname{\mathcal{S}}$, rather than prediction of the true response category probabilities where this method performs better than EN and GMR.

The average KL divergence and the average number of unique regression coefficient vectors estimated by GFMR for setting 1 when $\operatorname{\mathcal{L}}=\{(1,2),(2,3),(3,4)\}$ are reported in Table 2. These results show the same pattern that the results in Table 1 showed. While prediction accuracy decreased in this setting, we again see the GFMR using validation likelihood tuning is competitive with GMR. This decrease in prediction performance is expected because response categories are not naturally ordered in this data generating model. We also saw an increase in the average number of unique regression coefficient vectors estimated using GFMR with oracle tuning and GFMR with validation likelihood tuning when compared to the setting $\operatorname{\mathcal{L}}=\operatorname{\mathcal{S}}$.

To further investigate the number of unique regression coefficient vectors estimated by GFMR, we present a comparison between this quantity for the case when $\operatorname{\mathcal{L}}=\operatorname{\mathcal{S}}$ and $\operatorname{\mathcal{L}}=\{(1,2),(2,3),(3,4)\}$ in Table 3. We see that GFMR using validation likelihood performed poorly at selecting the true number of unique coefficient vectors estimated for both choices of $\operatorname{\mathcal{L}}$.

4.4 Results of Setting 2

In Table 4, we present the average KL divergence and the average number of unique regression coefficient vectors estimated by the GFMR method for the simulation using setting 2 when $\operatorname{\mathcal{L}}=\operatorname{\mathcal{S}}$. Similar to the patterns shown in setting 1 when $\operatorname{\mathcal{L}}=\operatorname{\mathcal{S}}$, on average the GFMR using validation likelihood tuning has a lower KL divergence than both EN and the GMR. For all values of $n$ and $\delta$, with the exception of $(\delta,n)=(0.1,50)$, both the oracle tuned GFMR and GFMR using validation likelihood tuning overselects the number of unique coefficient vectors. Again this may indicate the need for a different method for tuning parameter selection if our interest is model selection.

In Table 5, we present the average KL divergence and the average number of unique regression coefficients estimated by GFMR for the simulation using setting 2 when $\operatorname{\mathcal{L}}=\{(1,2),(2,3),(3,4)\}$. We see a similar result to the pattern observed Table 5, but the prediction accuracy here was better. This is expected because the response categories have a natural ordering in this data generating model.

To further investigate this comparison of unique regression coefficient vectors, Table 6 presents a comparison between the number of unique regression coefficient vectors estimated by the GFMR method in each of the 100 replications for the case when $\operatorname{\mathcal{L}}=\operatorname{\mathcal{S}}$ and the case when $\operatorname{\mathcal{L}}=\{(1,2),(2,3),(3,4)\}$. Just as in Section 4.3, these results show that GFMR using validation likelihood performs poorly at selecting the true number of unique coefficient vectors for both penalty sets, but again on average predicts the true response category probabilities better than EN and GMR.

5 Two-Step Method for Reducing Response Categories

The simulations presented in Section 4 show that our method using validation likelihood to select the tuning parameter performs well at predicting the true response category probabilities when compared to EN and GMR. The same simulations show that our method using validation likelihood does not perform well at model selection. To improve model selection performance, we propose an alternative approach for low dimensional settings. This is a two-step method that solves (4) for an increasing sequence of $\lambda$’s until there are two unique vectors in $(\hat{\beta}_{1},\ldots,\hat{\beta_{C}})$ . This sequence provides a set of candidate models with different response categories from which we will select the best by refitting using unpenalized maximum likelihood and computing the AIC. The selected model is the candidate model with the minimum AIC. We also compute the AIC for the edge case where all probabilities are equal to $1/C$.

We present a simulation to show the merits of this two-step approach to select the correct combined category model in multinomial logistic regression. The explanatory variables values are generated using the same procedure described in Section 4.2. In this simulation, $\beta^{*}_{1}=-\vec{\delta}$ and $\beta^{*}_{2}=\beta^{*}_{3}=\beta^{*}_{4}=\vec{0}$, and observed responses $y_{1},\ldots,y_{n}$ are a realization of $Y_{1},\ldots,Y_{n}$ defined by (1). We consider $(n,\delta)\in\{50,75\}\times\{1,3\}$.

To evaluate the proposed two-step method, we investigate its ability to detect the correct response category grouping. It is possible to detect the correct number of groups but the incorrect structure: we call this incorrect. It is also possible that the method selects more groups than it should, and if two of the groups were combined it would result in the correct structure: we call this one-step. Turning parameters were selected from a subset of $\{10^{-10},10^{-9.98},\ldots 10^{9.98},10^{10}\}$ and we set $\operatorname{\mathcal{L}}=\operatorname{\mathcal{S}}$. Since this simulation study has only 4 response categories, it was computationally feasible to compute the best AIC by searching over all possible response category combinations, which makes this exhaustive search a natural competitor.

In Table 7, we report the proportion of replications that the group structures of interest are selected for each $(n,\delta)$ combination for the two-step method. Table 8 presents the results for the exhaustive search. These tables show that for every of $n$ and $\delta$, the two-step method correctly picks the true response categories groups more than any other group structure. In particular, the proposed two-step method performs as well or better than exhaustive search for each $(n,\delta)$. These results show improvement in model selection performance when compared to the simulation results from Section 4, suggesting that if model selection is the interest the two step approach should be used.

6 Election Data Example

We analyze the dataset nes96 found in the CRAN package faraway (Faraway, 2014). The response variable, self-identified political affiliation of voters, has 7 levels: strong Democrat, weak Democrat, independent Democrat, independent, independent Republican, weak Republican, strong Republican. The explanatory variables are voter education level (categorical with 7 levels), voter income (categorical with 24 levels), and voter age (numerical). An investigation into both model selection based on validation likelihood tuning parameter selection and model selection based on using the two-step method was performed using $\operatorname{\mathcal{L}}=\operatorname{\mathcal{S}}$ and $\operatorname{\mathcal{L}}=\{(1,2),\ldots,(6,7)\}$.

Model selection for GFMR using validation likelihood tuning parameter selection was performed using 5-fold cross validation selecting the tuning parameter from the set

$\{10^{-10},10^{-9.98},\ldots,10^{9.98},10^{10}\}$. In Table 9, we present the response category combinations recommended by the GFMR regression coefficient estimates when $\operatorname{\mathcal{L}}=\operatorname{\mathcal{S}}$ and $\operatorname{\mathcal{L}}=\{(1,2),\ldots,(6,7)\}$. The results show for the case when $\operatorname{\mathcal{L}}=\operatorname{\mathcal{S}}$, GFMR does not combine any response categories. When $\operatorname{\mathcal{L}}=\{(1,2),\ldots,(6,7)\}$, the results show that independent Democrats, independent Republicans, and independents have the same estimated regression coefficient vectors.

We also show the results from the two-step approach proposed in Section 5. Both choices for $\operatorname{\mathcal{L}}$ resulted in a selected model with three response categories. These response category groups are shown in Table 10.