<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/more-flexible-ordinal-outcome-models/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://statcompute.wordpress.com/2018/08/28/more-flexible-ordinal-outcome-models/"> S+/R – Yet Another Blog in Statistical Computing</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>In the previous post (<a rel="nofollow" href="https://statcompute.wordpress.com/2018/08/26/adjacent-categories-and-continuation-ratio-logit-models-for-ordinal-outcomes" target="_blank">https://statcompute.wordpress.com/2018/08/26/adjacent-categories-and-continuation-ratio-logit-models-for-ordinal-outcomes</a>), we've shown alternative models for ordinal outcomes in addition to commonly used Cumulative Logit models under the proportional odds assumption, which are also known as Proportional Odds model. A potential drawback of Proportional Odds model is the lack of flexibility and the restricted assumption of proportional odds, of which the violation might lead to the model mis-specification. As a result, Cumulative Logit models with more flexible assumptions are called for. </p>
<p>In the example below, I will first show how to estimate credit ratings with a Cumulative Logit model under the proportional odds assumption with corporate financial performance measures, expressed as Logit(Y &lt;=j) = A_j – X * B, where A_j depends on the category j. </p>
<pre>
pkgs &lt;- list(&quot;maxLik&quot;, &quot;VGAM&quot;)
sapply(pkgs, require, character.only = T)
data(data_cr, envir = .GlobalEnv, package = &quot;mvord&quot;)
data_cr$RATING &lt;- pmax(data_cr$rater1, data_cr$rater2, na.rm = T)
x &lt;- c(&quot;LR&quot;, &quot;LEV&quot;, &quot;PR&quot;, &quot;RSIZE&quot;, &quot;BETA&quot;)
# LR   : LIQUIDITY RATIO
# LEV  : LEVERAGE RATIO
# PR   : PROFITABILITY RATIO
# RSIZE: LOG OF RELATIVE SIZE
# BETA : SYSTEMATIC RISK
y &lt;- &quot;RATING&quot;
df &lt;- data_cr[!is.na(data_cr[, y]), c(x, y)]
table(df[, y]) / length(df[, y])
#         A         B         C         D         E
# 0.1047198 0.1681416 0.3023599 0.2994100 0.1253687

### CUMULATIVE LOGIT MODEL ASSUMED PROPORTIONAL ODDS ###
# BELOW IS THE SIMPLER EQUIVALENT:
# vglm(RATING ~ LR + LEV + PR + RSIZE + BETA, data = df, family = cumulative(parallel = T))

ll1 &lt;- function(param) {
  plist &lt;- c(&quot;a_A&quot;, &quot;a_B&quot;, &quot;a_C&quot;, &quot;a_D&quot;, &quot;b_LR&quot;, &quot;b_LE&quot;, &quot;b_PR&quot;, &quot;b_RS&quot;, &quot;b_BE&quot;)
  sapply(1:length(plist), function(i) assign(plist[i], param[i], envir = .GlobalEnv))
  XB_A &lt;- with(df, a_A - (b_LR * LR + b_LE * LEV + b_PR * PR + b_RS * RSIZE + b_BE * BETA))
  XB_B &lt;- with(df, a_B - (b_LR * LR + b_LE * LEV + b_PR * PR + b_RS * RSIZE + b_BE * BETA))
  XB_C &lt;- with(df, a_C - (b_LR * LR + b_LE * LEV + b_PR * PR + b_RS * RSIZE + b_BE * BETA))
  XB_D &lt;- with(df, a_D - (b_LR * LR + b_LE * LEV + b_PR * PR + b_RS * RSIZE + b_BE * BETA))
  prob_A &lt;- 1 / (1 + exp(-XB_A))
  prob_B &lt;- 1 / (1 + exp(-XB_B)) - prob_A
  prob_C &lt;- 1 / (1 + exp(-XB_C)) - prob_A - prob_B
  prob_D &lt;- 1 / (1 + exp(-XB_D)) - prob_A - prob_B - prob_C
  prob_E &lt;- 1 - prob_A - prob_B - prob_C - prob_D
  CAT &lt;- data.frame(sapply(c(&quot;A&quot;, &quot;B&quot;, &quot;C&quot;, &quot;D&quot;, &quot;E&quot;), function(x) assign(x, df[, y] == x)))
  LH &lt;- with(CAT, (prob_A ^ A) * (prob_B ^ B) * (prob_C ^ C) * (prob_D ^ D) * (prob_E ^ E))
  return(sum(log(LH)))
}

start1 &lt;- c(a_A = 0, a_B = 2, a_C = 3, a_D = 4, b_LR = 0, b_LE = 0, b_PR = 0, b_RS = 0, b_BE = 0)
summary(m1 &lt;- maxLik(logLik = ll1, start = start1))
#     Estimate Std. error t value Pr(t)
#a_A  15.53765    0.77215  20.123  &lt;2e-16 ***
#a_B  18.26195    0.84043  21.729  &lt;2e-16 ***
#a_C  21.61729    0.94804  22.802  &lt;2e-16 ***
#a_D  25.88787    1.10522  23.423  &lt;2e-16 ***
#b_LR  0.29070    0.11657   2.494  0.0126 *
#b_LE  0.83977    0.07220  11.631  &lt;2e-16 ***
#b_PR -5.10955    0.35531 -14.381  &lt;2e-16 ***
#b_RS -2.18552    0.09982 -21.895  &lt;2e-16 ***
#b_BE  3.26811    0.21696  15.063  &lt;2e-16 ***
</pre>
<p>In the above output, the attribute &quot;liquidity ratioâ€� is somewhat less significant than the other, implying a potential opportunity for further improvements by relaxing the proportional odds assumption. As a result, I will try a different class of Cumulative Logit models, namely (unconstrained) Partial-Proportional Odds models, that would allow non-proportional odds for a subset of model attributes, e.g. LR in our case. Therefore, the formulation now becomes Logit(Y &lt;=j) = A_j – X * B – Z * G_j, where both A_j and G_j vary the category j. </p>
<pre>
### CUMULATIVE LOGIT MODEL ASSUMED UNCONSTRAINED PARTIAL-PROPORTIONAL ODDS ###
# BELOW IS THE SIMPLER EQUIVALENT:
# vglm(RATING ~ LR + LEV + PR + RSIZE + BETA, data = df, family = cumulative(parallel = F ~ LR))

ll2 &lt;- function(param) {
  plist &lt;- c(&quot;a_A&quot;, &quot;a_B&quot;, &quot;a_C&quot;, &quot;a_D&quot;, &quot;b_LRA&quot;, &quot;b_LRB&quot;, &quot;b_LRC&quot;, &quot;b_LRD&quot;, &quot;b_LE&quot;, &quot;b_PR&quot;, &quot;b_RS&quot;, &quot;b_BE&quot;)
  sapply(1:length(plist), function(i) assign(plist[i], param[i], envir = .GlobalEnv))
  XB_A &lt;- with(df, a_A - (b_LRA * LR + b_LE * LEV + b_PR * PR + b_RS * RSIZE + b_BE * BETA))
  XB_B &lt;- with(df, a_B - (b_LRB * LR + b_LE * LEV + b_PR * PR + b_RS * RSIZE + b_BE * BETA))
  XB_C &lt;- with(df, a_C - (b_LRC * LR + b_LE * LEV + b_PR * PR + b_RS * RSIZE + b_BE * BETA))
  XB_D &lt;- with(df, a_D - (b_LRD * LR + b_LE * LEV + b_PR * PR + b_RS * RSIZE + b_BE * BETA))
  prob_A &lt;- 1 / (1 + exp(-XB_A))
  prob_B &lt;- 1 / (1 + exp(-XB_B)) - prob_A
  prob_C &lt;- 1 / (1 + exp(-XB_C)) - prob_A - prob_B
  prob_D &lt;- 1 / (1 + exp(-XB_D)) - prob_A - prob_B - prob_C
  prob_E &lt;- 1 - prob_A - prob_B - prob_C - prob_D
  CAT &lt;- data.frame(sapply(c(&quot;A&quot;, &quot;B&quot;, &quot;C&quot;, &quot;D&quot;, &quot;E&quot;), function(x) assign(x, df[, y] == x)))
  LH &lt;- with(CAT, (prob_A ^ A) * (prob_B ^ B) * (prob_C ^ C) * (prob_D ^ D) * (prob_E ^ E))
  return(sum(log(LH)))
}

start2 &lt;- c(a_A = 0.1, a_B = 0.2, a_C = 0.3, a_D = 0.4, b_LRA = 0, b_LRB = 0, b_LRC = 0, b_LRD = 0, b_LE = 0, b_PR = 0, b_RS = 0, b_BE = 0)
summary(m2 &lt;- maxLik(logLik = ll2, start = start2))
#Estimates:
#      Estimate Std. error t value Pr(t)
#a_A   15.30082    0.83936  18.229  &lt;2e-16 ***
#a_B   18.14795    0.81325  22.315  &lt;2e-16 ***
#a_C   21.72469    0.89956  24.150  &lt;2e-16 ***
#a_D   25.92697    1.07749  24.062  &lt;2e-16 ***
#b_LRA  0.12442    0.30978   0.402  0.6880
#b_LRB  0.21127    0.20762   1.018  0.3089
#b_LRC  0.36097    0.16687   2.163  0.0305 *
#b_LRD  0.31404    0.22090   1.422  0.1551
#b_LE   0.83949    0.07155  11.733  &lt;2e-16 ***
#b_PR  -5.09891    0.35249 -14.465  &lt;2e-16 ***
#b_RS  -2.18589    0.09540 -22.913  &lt;2e-16 ***
#b_BE   3.26529    0.20993  15.554  &lt;2e-16 ***
</pre>
<p>As shown above, under the partial-proportional odds assumption, there are 4 parameters estimated for LR, three of which are not significant and therefore the additional flexibility is not justified. In fact, AIC of the 2nd model (AIC = 1103.60) is even higher than AIC of the 1st model (AIC = 1098.18). </p>
<p>In light of the above observation, I will introduce the 3rd model, which is known as the Constrained Partial-Proportional Odds model and expressed as Logit(Y &lt;=j) = A_j – X * B – Z * G * gamma_j, where A_j and gamma_j vary the category j. It is worth pointing out that gamma_j is a pre-specified fixed scalar and does not need to be estimated. Based on the unconstrained model outcome, we can set gamma_1 = 1, gamma_2 = 2, and gamma_3 = gamma_4 = 3 for LR in our case. </p>
<pre>
### CUMULATIVE LOGIT MODEL ASSUMED CONSTRAINED PARTIAL-PROPORTIONAL ODDS ###

ll3 &lt;- function(param) {
  plist &lt;- c(&quot;a_A&quot;, &quot;a_B&quot;, &quot;a_C&quot;, &quot;a_D&quot;, &quot;b_LR&quot;, &quot;b_LE&quot;, &quot;b_PR&quot;, &quot;b_RS&quot;, &quot;b_BE&quot;)
  sapply(1:length(plist), function(i) assign(plist[i], param[i], envir = .GlobalEnv))
  gamma &lt;- c(1, 2, 3, 3)
  XB_A &lt;- with(df, a_A - (gamma[1] * b_LR * LR + b_LE * LEV + b_PR * PR + b_RS * RSIZE + b_BE * BETA))
  XB_B &lt;- with(df, a_B - (gamma[2] * b_LR * LR + b_LE * LEV + b_PR * PR + b_RS * RSIZE + b_BE * BETA))
  XB_C &lt;- with(df, a_C - (gamma[3] * b_LR * LR + b_LE * LEV + b_PR * PR + b_RS * RSIZE + b_BE * BETA))
  XB_D &lt;- with(df, a_D - (gamma[4] * b_LR * LR + b_LE * LEV + b_PR * PR + b_RS * RSIZE + b_BE * BETA))
  prob_A &lt;- 1 / (1 + exp(-XB_A))
  prob_B &lt;- 1 / (1 + exp(-XB_B)) - prob_A
  prob_C &lt;- 1 / (1 + exp(-XB_C)) - prob_A - prob_B
  prob_D &lt;- 1 / (1 + exp(-XB_D)) - prob_A - prob_B - prob_C
  prob_E &lt;- 1 - prob_A - prob_B - prob_C - prob_D
  CAT &lt;- data.frame(sapply(c(&quot;A&quot;, &quot;B&quot;, &quot;C&quot;, &quot;D&quot;, &quot;E&quot;), function(x) assign(x, df[, y] == x)))
  LH &lt;- with(CAT, (prob_A ^ A) * (prob_B ^ B) * (prob_C ^ C) * (prob_D ^ D) * (prob_E ^ E))
  return(sum(log(LH)))
}

start3 &lt;- c(a_A = 1, a_B = 2, a_C = 3, a_D = 4, b_LR = 0.1, b_LE = 0, b_PR = 0, b_RS = 0, b_BE = 0)
summary(m3 &lt;- maxLik(logLik = ll3, start = start3))
#Estimates:
#     Estimate Std. error t value Pr(t)
#a_A  15.29442    0.60659  25.214 &lt; 2e-16 ***
#a_B  18.18220    0.65734  27.660 &lt; 2e-16 ***
#a_C  21.70599    0.75181  28.872 &lt; 2e-16 ***
#a_D  25.98491    0.88104  29.493 &lt; 2e-16 ***
#b_LR  0.11351    0.04302   2.638 0.00833 **
#b_LE  0.84012    0.06939  12.107 &lt; 2e-16 ***
#b_PR -5.10025    0.33481 -15.233 &lt; 2e-16 ***
#b_RS -2.18708    0.08118 -26.940 &lt; 2e-16 ***
#b_BE  3.26689    0.19958  16.369 &lt; 2e-16 ***
</pre>
<p>As shown above, after the introduction of gamma_j as the constrained scalar, the statistical significance of LR has been improved with a slightly lower AIC = 1097.64. </p>
<p>To be complete, I'd like to mention the last model today, which is named the Stereotype model. The idea of Stereotype models is very similar to the idea of adjacent-categories models and is to estimate Log(Y = j / Y = j+1) or more often Log(Y = j / Y = j_c), where C represents a baseline category. However, the right-hand side is expressed as Log(…) = A_j – (X * B) * phi_j, where phi_j is a hyper-parameter such that phi_1 = 1 > phi_2…> phi_max = 0. As a result, the coefficient of each model attribute could also vary by the category j, introducing more flexibility at the cost of being difficult to estimate. </p>
<pre>
### STEREOTYPE MODEL ###
# BELOW IS THE SIMPLER EQUIVALENT:
# rrvglm(sapply(c(&quot;A&quot;, &quot;B&quot;, &quot;C&quot;, &quot;D&quot;, &quot;E&quot;), function(x) df[, y] == x)~ LR + LEV + PR + RSIZE + BETA, multinomial, data = df)

ll4 &lt;- function(param) {
  plist &lt;- c(&quot;a_A&quot;, &quot;a_B&quot;, &quot;a_C&quot;, &quot;a_D&quot;, &quot;b_LR&quot;, &quot;b_LE&quot;, &quot;b_PR&quot;, &quot;b_RS&quot;, &quot;b_BE&quot;, &quot;phi_B&quot;, &quot;phi_C&quot;, &quot;phi_D&quot;)
  sapply(1:length(plist), function(i) assign(plist[i], param[i], envir = .GlobalEnv))
  XB_A &lt;- with(df, a_A - (b_LR * LR + b_LE * LEV + b_PR * PR + b_RS * RSIZE + b_BE * BETA))
  XB_B &lt;- with(df, a_B - phi_B * (b_LR * LR + b_LE * LEV + b_PR * PR + b_RS * RSIZE + b_BE * BETA))
  XB_C &lt;- with(df, a_C - phi_C * (b_LR * LR + b_LE * LEV + b_PR * PR + b_RS * RSIZE + b_BE * BETA))
  XB_D &lt;- with(df, a_D - phi_D * (b_LR * LR + b_LE * LEV + b_PR * PR + b_RS * RSIZE + b_BE * BETA))
  prob_A &lt;- exp(XB_A) / (exp(XB_A) + exp(XB_B) + exp(XB_C) + exp(XB_D) + 1)
  prob_B &lt;- exp(XB_B) / (exp(XB_A) + exp(XB_B) + exp(XB_C) + exp(XB_D) + 1)
  prob_C &lt;- exp(XB_C) / (exp(XB_A) + exp(XB_B) + exp(XB_C) + exp(XB_D) + 1)
  prob_D &lt;- exp(XB_D) / (exp(XB_A) + exp(XB_B) + exp(XB_C) + exp(XB_D) + 1)
  prob_E &lt;- 1 - prob_A - prob_B - prob_C - prob_D
  CAT &lt;- data.frame(sapply(c(&quot;A&quot;, &quot;B&quot;, &quot;C&quot;, &quot;D&quot;, &quot;E&quot;), function(x) assign(x, df[, y] == x)))
  LH &lt;- with(CAT, (prob_A ^ A) * (prob_B ^ B) * (prob_C ^ C) * (prob_D ^ D) * (prob_E ^ E))
  return(sum(log(LH)))
}

start4 &lt;- c(a_A = 1, a_B = 2, a_C = 3, a_D = 4, b_LR = 0.1, b_LE = 0, b_PR = 0, b_RS = 0, b_BE = 0, phi_B = 0.1, phi_C = 0.2, phi_D = 0.3)
summary(m4 &lt;- maxLik(logLik = ll4, start = start4))
#Estimates:
#       Estimate Std. error t value Pr(t)
#a_A    67.73429    2.37424  28.529  &lt;2e-16 ***
#a_B    55.86469    1.94442  28.731  &lt;2e-16 ***
#a_C    41.27477    1.47960  27.896  &lt;2e-16 ***
#a_D    22.24244    1.83137  12.145  &lt;2e-16 ***
#b_LR    0.86975    0.37481   2.320  0.0203 *
#b_LE    2.79215    0.23373  11.946  &lt;2e-16 ***
#b_PR  -16.66836    1.17569 -14.178  &lt;2e-16 ***
#b_RS   -7.24921    0.33460 -21.665  &lt;2e-16 ***
#b_BE   10.57411    0.72796  14.526  &lt;2e-16 ***
#phi_B   0.77172    0.03155  24.461  &lt;2e-16 ***
#phi_C   0.52806    0.03187  16.568  &lt;2e-16 ***
#phi_D   0.26040    0.02889   9.013  &lt;2e-16 ***
</pre>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://statcompute.wordpress.com/2018/08/28/more-flexible-ordinal-outcome-models/"> S+/R – Yet Another Blog in Statistical Computing</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/wzHUzXEHM7I" height="1"></a></div>