Models for Polytomous Responses
Shreejit
6/5/2020
#install.packages("car")
#install.packages("VGAM")For a target with Ordinal Values
Proportional Odds Model
Unless the response variable has numeric values, it is important to ensure that it has been defined as an ordered factor (using ordered()). In the Arthritis data, the response, Improved was setup this way, as we can check by printing some of the values
library(vcd)## Loading required package: gridlibrary(grid)
data("Arthritis", package="vcd")
head(Arthritis$Improved, 8)## [1] Some None None Marked Marked Marked None Marked
## Levels: None < Some < Markedlibrary(MASS)
arth.polr <- polr(Improved ~ Sex + Treatment + Age,data=Arthritis, Hess=TRUE)
summary(arth.polr)## Call:
## polr(formula = Improved ~ Sex + Treatment + Age, data = Arthritis,
## Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## SexMale -1.25168 0.54636 -2.291
## TreatmentTreated 1.74529 0.47589 3.667
## Age 0.03816 0.01842 2.072
##
## Intercepts:
## Value Std. Error t value
## None|Some 2.5319 1.0571 2.3952
## Some|Marked 3.4309 1.0912 3.1443
##
## Residual Deviance: 145.4579
## AIC: 155.4579Summary gives us the coefficients, Intercepts labeled by the cutpoint on the ordinal response but no p-values. The car::Anova() method gives the appropriate tests
library(car)## Loading required package: carDataAnova(arth.polr)## Analysis of Deviance Table (Type II tests)
##
## Response: Improved
## LR Chisq Df Pr(>Chisq)
## Sex 5.6880 1 0.0170812 *
## Treatment 14.7095 1 0.0001254 ***
## Age 4.5715 1 0.0325081 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In R, the PO and NPO models can be readily contrasted by fitting them both using vglm() in the VGAM package. This defines the cumulative family of models and allows a parallel option. With parallel=TRUE, this is equivalent to the polr() model, except that the signs of the coefficients are reversed
#install.packages("VGAM")
library(VGAM)## Loading required package: stats4## Loading required package: splines##
## Attaching package: 'VGAM'## The following object is masked from 'package:car':
##
## logitarth.po <- vglm(Improved ~ Sex + Treatment + Age, data=Arthritis,family = cumulative(parallel=TRUE))
arth.po##
## Call:
## vglm(formula = Improved ~ Sex + Treatment + Age, family = cumulative(parallel = TRUE),
## data = Arthritis)
##
##
## Coefficients:
## (Intercept):1 (Intercept):2 SexMale TreatmentTreated
## 2.53199006 3.43098785 1.25167130 -1.74530408
## Age
## -0.03816279
##
## Degrees of Freedom: 168 Total; 163 Residual
## Residual deviance: 145.4579
## Log-likelihood: -72.72896The more general NPO model can be fit using parallel=FALSE.
arth.npo <- vglm(Improved ~ Sex + Treatment + Age, data=Arthritis, family = cumulative(parallel=FALSE))
arth.npo##
## Call:
## vglm(formula = Improved ~ Sex + Treatment + Age, family = cumulative(parallel = FALSE),
## data = Arthritis)
##
##
## Coefficients:
## (Intercept):1 (Intercept):2 SexMale:1
## 2.61853875 3.43117450 1.50982704
## SexMale:2 TreatmentTreated:1 TreatmentTreated:2
## 0.86643395 -1.83692928 -1.70401146
## Age:1 Age:2
## -0.04086648 -0.03729421
##
## Degrees of Freedom: 168 Total; 160 Residual
## Residual deviance: 143.5741
## Log-likelihood: -71.78703The VGAM package defines a coef() method that can print the coefficients in a more readable matrix form giving the category cutpoints
coef(arth.po, matrix=TRUE)## logitlink(P[Y<=1]) logitlink(P[Y<=2])
## (Intercept) 2.53199006 3.43098785
## SexMale 1.25167130 1.25167130
## TreatmentTreated -1.74530408 -1.74530408
## Age -0.03816279 -0.03816279coef(arth.npo, matrix=TRUE)## logitlink(P[Y<=1]) logitlink(P[Y<=2])
## (Intercept) 2.61853875 3.43117450
## SexMale 1.50982704 0.86643395
## TreatmentTreated -1.83692928 -1.70401146
## Age -0.04086648 -0.03729421In most cases, nested models can be tested using an anova() method, but the VGAM pack- age has not implemented this for “vglm” objects. Instead, it provides an analogous function, lrtest()
VGAM::lrtest(arth.npo, arth.po)## Likelihood ratio test
##
## Model 1: Improved ~ Sex + Treatment + Age
## Model 2: Improved ~ Sex + Treatment + Age
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 160 -71.787
## 2 163 -72.729 3 1.8838 0.5969The LR test can be also calculated as “manually” shown below using the difference in residual deviance for the two models
tab <- cbind(Deviance = c(deviance(arth.npo), deviance(arth.po)),df = c(df.residual(arth.npo), df.residual(arth.po)))
tab <- rbind(tab, diff(tab))
rownames(tab) <- c("GenLogit", "PropOdds", "LR test")
tab <- cbind(tab, pvalue=1-pchisq(tab[,1], tab[,2]))
tab## Deviance df pvalue
## GenLogit 143.574067 160 0.8196626
## PropOdds 145.457915 163 0.8343527
## LR test 1.883849 3 0.5968607The vglm() can also fit partial proportional odds models, by specifying a formula giving the terms for which the PO assumption should be taken as TRUE or FALSE. Here we illustrate this using parallel=FALSE ~ Sex,tofitseparateslopesformalesandfemales,butparallellinesforthe otherpredictors.Thesamemodelwouldbefitusingparallel=TRUE ~ Treatment + Age.
arth.ppo <- vglm(Improved ~ Sex + Treatment + Age, data=Arthritis,family = cumulative(parallel=FALSE ~ Sex))
coef(arth.ppo, matrix=TRUE)## logitlink(P[Y<=1]) logitlink(P[Y<=2])
## (Intercept) 2.54245216 3.61556073
## SexMale 1.48333554 0.86736192
## TreatmentTreated -1.77574178 -1.77574178
## Age -0.03962233 -0.03962233Graphical Assesment of Proportional Odds
#install.packages("polspline",dependencies = TRUE)
#install.packages("rms",dependencies = TRUE)
library("rms")## Loading required package: Hmisc## Loading required package: lattice## Loading required package: survival## Loading required package: Formula## Loading required package: ggplot2##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:base':
##
## format.pval, units## Loading required package: SparseM##
## Attaching package: 'SparseM'## The following object is masked from 'package:base':
##
## backsolve##
## Attaching package: 'rms'## The following objects are masked from 'package:VGAM':
##
## calibrate, lrtest## The following objects are masked from 'package:car':
##
## Predict, vifarth.po2 <- lrm(Improved ~ Sex + Treatment + Age, data=Arthritis)
arth.po2## Logistic Regression Model
##
## lrm(formula = Improved ~ Sex + Treatment + Age, data = Arthritis)
##
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 84 LR chi2 24.46 R2 0.291 C 0.750
## None 42 d.f. 3 g 1.335 Dxy 0.500
## Some 14 Pr(> chi2) <0.0001 gr 3.801 gamma 0.503
## Marked 28 gp 0.280 tau-a 0.309
## max |deriv| 1e-07 Brier 0.187
##
## Coef S.E. Wald Z Pr(>|Z|)
## y>=Some -2.5320 1.0570 -2.40 0.0166
## y>=Marked -3.4310 1.0911 -3.14 0.0017
## Sex=Male -1.2517 0.5464 -2.29 0.0220
## Treatment=Treated 1.7453 0.4759 3.67 0.0002
## Age 0.0382 0.0184 2.07 0.0382
## op <- par(mfrow=c(1,3))
plot.xmean.ordinaly(Improved ~ Sex + Treatment + Age, data=Arthritis,lwd=2, pch=16, subn=FALSE)
par(op)Visualizing results for Proportional Odds Model
arth.fitp <- cbind(Arthritis,predict(arth.polr, type="probs"))
head(arth.fitp) ## ID Treatment Sex Age Improved None Some Marked
## 1 57 Treated Male 27 Some 0.7326185 0.1380586 0.1293229
## 2 46 Treated Male 29 None 0.7174048 0.1444325 0.1381627
## 3 77 Treated Male 30 None 0.7096042 0.1476259 0.1427699
## 4 17 Treated Male 32 Marked 0.6936286 0.1540035 0.1523679
## 5 36 Treated Male 46 Marked 0.5702499 0.1950359 0.2347142
## 6 23 Treated Male 58 Marked 0.4563432 0.2171302 0.3265266library(reshape2)
plotdat <- melt(arth.fitp,id.vars = c("Sex", "Treatment", "Age", "Improved"),measure.vars=c("None", "Some", "Marked"),
variable.name = "Level",
value.name = "Probability")
head(plotdat)## Sex Treatment Age Improved Level Probability
## 1 Male Treated 27 Some None 0.7326185
## 2 Male Treated 29 None None 0.7174048
## 3 Male Treated 30 None None 0.7096042
## 4 Male Treated 32 Marked None 0.6936286
## 5 Male Treated 46 Marked None 0.5702499
## 6 Male Treated 58 Marked None 0.4563432library(ggplot2)
library(directlabels)
gg <- ggplot(plotdat, aes(x = Age, y = Probability, colour = Level)) +
geom_line(size=2.5) + theme_bw() + xlim(10,80) +
geom_point(color="black", size=1.5) +
facet_grid(Sex ~ Treatment) #,labeller = function(x, y) sprintf("%s = %s", x, y)
direct.label(gg)
Although we now have three response curves in each panel, this plot is relatively easy to un- derstand: (a) In each panel, the probability of no improvement decreases with age, while that for marked improvement increases. (b) It is easy to compare the placebo and treated groups in each row, showing that no improvement decreases, while marked improvement increases with the active treatment. (On the other hand, this layout makes it harder to compare panels vertically for males and females in each condition.) (c) The points show where the observations are located in each panel; so, we can see that the data is quite thin for males given the placebo
Effect Plots
#install.packages("effects")
library(effects)## Use the command
## lattice::trellis.par.set(effectsTheme())
## to customize lattice options for effects plots.
## See ?effectsTheme for details.plot(Effect("Age", arth.polr))
plot(Effect("Age", arth.polr), style='stacked',key.args=list(x=.55, y=.9))
plot(Effect(c("Treatment", "Sex", "Age"), arth.polr),style="stacked", key.arg=list(x=.8, y=.9))
Latent variable effect plot for the effects of Treatment and Age in the Arthritis data
plot(Effect(c("Treatment", "Age"), arth.polr, latent=TRUE), lwd=3)
EXAMPLE 8.1: Women’s labor force participation
library(car)
rm(Womenlf)
data("Womenlf", package="car")
#head(Womenlf)Womenlf <- WomenlfIn this example, it makes sense to consider a first dichotomy (working) between women who are not working, vs. those who are (full time or part time). A second dichotomy (fulltime) contrasts full time work vs. part time work, among those women who are working at least part time. These two binary variables are created in the data frame using the recode() function from the car package
#Create dichotomies// using ifelse. Recode is another option which was not working in this case
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:Hmisc':
##
## src, summarize## The following object is masked from 'package:car':
##
## recode## The following object is masked from 'package:MASS':
##
## select## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionworking = ifelse(Womenlf$partic=='not.work',"no","yes")
fulltime = ifelse(Womenlf$partic=='fulltime',"yes",
ifelse(Womenlf$partic=='parttime',"no",NA)
)
Womenlf = cbind(Womenlf,working,fulltime)
#Womenlf = cbind(Womenlf,fulltime)
head(Womenlf)## partic hincome children region working fulltime
## 1 not.work 15 present Ontario no <NA>
## 2 not.work 13 present Ontario no <NA>
## 3 not.work 45 present Ontario no <NA>
## 4 not.work 23 present Ontario no <NA>
## 5 not.work 19 present Ontario no <NA>
## 6 not.work 7 present Ontario no <NA>with(Womenlf, table(partic, working))## working
## partic no yes
## fulltime 0 66
## not.work 155 0
## parttime 0 42with(Womenlf, table(partic, fulltime, useNA="ifany"))## fulltime
## partic no yes <NA>
## fulltime 0 66 0
## not.work 0 0 155
## parttime 42 0 0We proceed to fit two separate binary logistic regression models for the derived dichotomous variables. For the working dichotomy, we get the following results:
mod.working <- glm(working ~ hincome + children, family=binomial,data=Womenlf)
summary(mod.working)##
## Call:
## glm(formula = working ~ hincome + children, family = binomial,
## data = Womenlf)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6767 -0.8652 -0.7768 0.9292 1.9970
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.33583 0.38376 3.481 0.0005 ***
## hincome -0.04231 0.01978 -2.139 0.0324 *
## childrenpresent -1.57565 0.29226 -5.391 7e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 356.15 on 262 degrees of freedom
## Residual deviance: 319.73 on 260 degrees of freedom
## AIC: 325.73
##
## Number of Fisher Scoring iterations: 4mod.fulltime <- glm(fulltime ~ hincome + children, family=binomial,data=Womenlf)
summary(mod.fulltime)##
## Call:
## glm(formula = fulltime ~ hincome + children, family = binomial,
## data = Womenlf)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4047 -0.8678 0.3949 0.6213 1.7641
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.47777 0.76711 4.534 5.80e-06 ***
## hincome -0.10727 0.03915 -2.740 0.00615 **
## childrenpresent -2.65146 0.54108 -4.900 9.57e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 144.34 on 107 degrees of freedom
## Residual deviance: 104.49 on 105 degrees of freedom
## (155 observations deleted due to missingness)
## AIC: 110.49
##
## Number of Fisher Scoring iterations: 5Although these were fit separately, we can view this as a combined model for the three-level response, with the following coefficients
cbind(working=coef(mod.working), fulltime=coef(mod.fulltime))## working fulltime
## (Intercept) 1.33582979 3.4777735
## hincome -0.04230843 -0.1072679
## childrenpresent -1.57564843 -2.6514557For both dichotomies, increasing income of the husband and the presence of young children decrease the log odds of a greater level of work. However, for those women who are working the effects of husband’s income and and children are greater on the choice between full time and part time work than they are for all women on the choice between working and not workin
LRtest <- function(model) c(LRchisq=(model$null.deviance - model$deviance),df=(model$df.null - model$df.residual))
tab <- rbind(working=LRtest(mod.working),fulltime=LRtest(mod.fulltime))
tab <- rbind(tab, All = colSums(tab))
tab <- cbind(tab, pvalue = 1- pchisq(tab[,1], tab[,2]))
tab## LRchisq df pvalue
## working 36.41835 2 1.235536e-08
## fulltime 39.84682 2 2.225215e-09
## All 76.26518 4 1.110223e-15Anova(mod.working)## Analysis of Deviance Table (Type II tests)
##
## Response: working
## LR Chisq Df Pr(>Chisq)
## hincome 4.8264 1 0.02803 *
## children 31.3229 1 2.185e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(mod.fulltime)## Analysis of Deviance Table (Type II tests)
##
## Response: fulltime
## LR Chisq Df Pr(>Chisq)
## hincome 8.981 1 0.002728 **
## children 32.136 1 1.437e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1predictors <- expand.grid(hincome=1:50,children=c('absent', 'present'))
fit <- data.frame(predictors,
p.working = predict(mod.working, predictors, type='response'),
p.fulltime = predict(mod.fulltime, predictors, type='response'),
l.working = predict(mod.working, predictors, type='link'),
l.fulltime = predict(mod.fulltime, predictors, type='link')
)
print(some(fit, 5), digits=3)## hincome children p.working p.fulltime l.working l.fulltime
## 18 18 absent 0.640 0.8245 0.574 1.547
## 29 29 absent 0.527 0.5907 0.109 0.367
## 34 34 absent 0.474 0.4578 -0.103 -0.169
## 37 37 absent 0.443 0.3796 -0.230 -0.491
## 96 46 present 0.101 0.0162 -2.186 -4.108fit <- within(fit, { full <- p.working * p.fulltime
part <- p.working * (1 - p.fulltime)
not <- 1 - p.working
})op <- par(mfrow=c(1,2), mar=c(5,4,4,1)+.1)
Hinc <- 1:max(fit$hincome)
for ( kids in c("absent", "present") )
{
dat <- subset(fit,children==kids)
plot(range(Hinc),c(0,1),type="n",cex.lab=1.25,
xlab="Husbands Income",ylab="Fitted Probability",
main=paste("Children",kids))
lines(Hinc, dat$not, lwd=3, col="black", lty=1)
lines(Hinc, dat$part, lwd=3, col="blue", lty=2)
lines(Hinc, dat$full, lwd=3, col="red", lty=5)
if(kids=="absent")
{
legend("topright", lty=c(1,2,5), lwd=3, col=c("black", "blue", "red"),
legend=c('not working', 'part-time', 'full-time'))
}
} 
par(op)We can see how that the decision not to work outside the home increases strongly with husband’s income, and is higher when there are children present. As well, among working women, the decision to work full time as opposed to part time decreases strongly with husband’s income, and is less likely with young children. Similarly, we plot the fitted logits for the two dichotomies in l.working and l.fulltime as shown below:
op <- par(mfrow=c(1,2), mar=c(5,4,1,1)+.1)
for ( kids in c("absent", "present") )
{
dat <- subset(fit,children==kids)
plot(range(Hinc),c(0,1),type="n",cex.lab=1.25,
xlab="Husbands Income",ylab="Fitted Log Odds")
lines(Hinc, dat$l.working, lwd=3, col="black", lty=1)
lines(Hinc, dat$l.fulltime, lwd=3, col="blue", lty=2)
text(25, 4.5, paste("Children", kids), cex=1.4)
if (kids=="absent")
{
legend("bottomleft", lty=1:2, lwd=3, col=c("black", "blue"),
legend=c('working', 'full-time'))
}
}
These Graphs didn’t come out as expected
Generalized Logit Model
Response : Women working full time, part time or not working Variables: Husbands Income (Continuous), Children Present vs. Absent (Categorical)
levels(Womenlf$partic)## [1] "fulltime" "not.work" "parttime"# choose not working as baseline category
Womenlf$partic <- relevel(Womenlf$partic, ref="not.work")We fit the main effects model for husband’s income and children as follows. As we did with polr() (Section 8.1), specifying Hess=TRUE saves the Hessian and facilitates calculation of standard errors and hypothesis tests
library(nnet)
wlf.multinom <- multinom(partic ~ hincome + children,data=Womenlf, Hess=TRUE)## # weights: 12 (6 variable)
## initial value 288.935032
## iter 10 value 211.454772
## final value 211.440963
## convergedThe summary() method for “multinom” objects doesn’t calculate test statistics for the esti- mated coefficients by default. The option Wald=TRUE produces Wald z-test statistics, calculated as z = β/SE(β)
summary(wlf.multinom, Wald=TRUE)## Call:
## multinom(formula = partic ~ hincome + children, data = Womenlf,
## Hess = TRUE)
##
## Coefficients:
## (Intercept) hincome childrenpresent
## fulltime 1.982842 -0.097232073 -2.55860537
## parttime -1.432321 0.006893838 0.02145558
##
## Std. Errors:
## (Intercept) hincome childrenpresent
## fulltime 0.4841789 0.02809599 0.3621999
## parttime 0.5924627 0.02345484 0.4690352
##
## Value/SE (Wald statistics):
## (Intercept) hincome childrenpresent
## fulltime 4.095266 -3.4607098 -7.06407045
## parttime -2.417573 0.2939197 0.04574407
##
## Residual Deviance: 422.8819
## AIC: 434.8819Notice that the coefficients, their standard errors and the Wald test z values are printed in separate tables. The first line in each table pertains to the logit comparing full time work with the not working reference level; the second line compares part time work against not working
For those who like p-values for significance tests, you can calculate these from the results re- turned by the summary() method in the Wald.ratios component, using the standard normal asymptotic approximation
stats <- summary(wlf.multinom, Wald=TRUE)
z <- stats$Wald.ratios
p <- 2 * (1 - pnorm(abs(z)))
zapsmall(p)## (Intercept) hincome childrenpresent
## fulltime 0.0000422 0.0005388 0.0000000
## parttime 0.0156244 0.7688193 0.9635142wlf.multinom2 <- multinom(partic ~ hincome * children, data=Womenlf, Hess=TRUE)## # weights: 15 (8 variable)
## initial value 288.935032
## iter 10 value 210.797079
## final value 210.714841
## convergedAnova(wlf.multinom2)## Analysis of Deviance Table (Type II tests)
##
## Response: partic
## LR Chisq Df Pr(>Chisq)
## hincome 15.153 2 0.0005123 ***
## children 63.559 2 1.579e-14 ***
## hincome:children 1.452 2 0.4837815
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Full model plots of the fitted values can be plotted as shown earlier in Example 8.1: obtain the fitted values over a grid of the predictors and plot these
predictors <- expand.grid(hincome=1:50,
children=c('absent', 'present'))
fit <- data.frame(predictors,predict(wlf.multinom, predictors, type='probs'))Plotting these fitted values gives the plot shown in Figure 8.12.
op <- par(mfrow=c(1,2), mar=c(5,4,4,1)+.1)
Hinc <- 1:max(fit$hincome)
for ( kids in c("absent", "present") ) {
dat <- subset(fit, children==kids)
plot( range(Hinc), c(0,1), type="n", cex.lab=1.25,
xlab="Husband's Income", ylab='Fitted Probability',
main = paste("Children", kids))
lines(Hinc, dat$not.work, lwd=3, col="black", lty=1)
lines(Hinc, dat$parttime, lwd=3, col="blue", lty=2)
lines(Hinc, dat$fulltime, lwd=3, col="red", lty=5)
if (kids=="absent") {
legend("topright", lty=c(1,2,5), lwd=3, col=c("black", "blue", "red"),legend=c('not working', 'part-time','full-time'))
}
} 
The effects package has special methods for “multinom” models. It treats the response levels in the order given by levels(), so before plotting we use ordered() to arrange levels in their nat- ural order. The update() method provides a simple way to get a new fitted model; in the call, the modelformula. ~ .meanstofitthesamemodelasbefore,i.e.,partic ~ hincome + children
Womenlf$partic <- ordered(Womenlf$partic,
levels=c('not.work', 'parttime', 'fulltime'))
wlf.multinom <- update(wlf.multinom, . ~ .)## # weights: 12 (6 variable)
## initial value 288.935032
## iter 10 value 211.454772
## final value 211.440963
## convergedAsillustratedearlier,youcanuseplot(allEffects(model), …)toplotallthehigh- order terms in the model, either with separate curves for each response level (style=“lines”) or as cumulative filled polygons (style=“stacked”). Here, we simply plot the effects for the combinations of husband’s income and children in stacked style, giving a plot (Figure 8.13) that is analogous to the full-model plot shown in Figure 8.12.
plot(Effect(c("hincome", "children"), wlf.multinom),
style="stacked", key.args=list(x=.05, y=.9))