无序多分类Logistic回归

• 定义：研究的因变量类型为多分类无序资料，应采用多分类无序Logistic回归模型分析。（有序多分类Logistic回归模型在调整连接函数后仍不满足平行线假设，可考虑采用多分类无序Logistic回归模型）。
  无序多分类logistic回归分析，是用于研究自变量 X 与 因变量 Y（分类变量）的一种多元回归分析法。自变量X 可以为分类变量或者连续型变量，因变量Y 为无序多分类变量，即不考虑因变量之间的等级信息，例如：职业可分为医生、教师、工程师、工人等。
  多分类无序Logistic回归模型，以某一水平因变量为参照水平，与其他水平相比建立（水平数-1）个广义Logistic模型

1. 模型构建

library(foreign)
library(nnet)
library(ggplot2)
library(reshape2)
ml <- read.dta("hsbdemo.dta") #hsbdemo.dta为示例数据
with(ml, table(ses, prog))
with(ml, do.call(rbind, tapply(write, prog, function(x) c(M = mean(x), SD = sd(x)))))
ml$prog2 <- relevel(ml$prog, ref = "academic")
test <- multinom(prog2 ~ ses + write, data = ml)
summary(test)

2. 统计检验

# 2-tailed z test
z <- summary(test)$coefficients/summary(test)$standard.errors 
z
# (Intercept)  sesmiddle   seshigh     write
# general     2.445214 -1.2018081 -2.261334 -2.705562
# vocation    4.484769  0.6116747 -1.649967 -5.112689
p <- (1 - pnorm(abs(z), 0, 1)) * 2

3. 模型预测

# extract the coefficients from the model and exponentiate
exp(coef(test))
# (Intercept) sesmiddle   seshigh     write
# general     17.32582 0.5866769 0.3126026 0.9437172
# vocation   184.61262 1.3382809 0.3743123 0.8926116
head(pp <- fitted(test))

dses <- data.frame(ses = c("low", "middle", "high"), write = mean(ml$write))
predict(test, newdata = dses, "probs")

dwrite <- data.frame(ses = rep(c("low", "middle", "high"), each = 41), write = rep(c(30:70),3))

# store the predicted probabilities for each value of ses and write
pp.write <- cbind(dwrite, predict(test, newdata = dwrite, type = "probs", se = TRUE))

# calculate the mean probabilities within each level of ses
by(pp.write[, 3:5], pp.write$ses, colMeans)

#melt data set to long for ggplot2
lpp <- melt(pp.write, id.vars = c("ses", "write"), value.name = "probability")
head(lpp)  # view first few rows

ggplot(lpp, aes(x = write, y = probability, colour = ses)) + geom_line() + facet_grid(variable ~., scales = "free")

等级（有序多分类）Logistic回归

• 有序多分类的Logistic回归原理是将因变量的多个分类依次分割为多个二元的Logistic回归，例如本例中因变量“本国的税收过高”的赞同程度有4个等级，分析时拆分为三个二元Logistic回归，分别为(0 vs 1+2+3) 、(0+1 vs 2+3)、(0+1+2 vs 3)，均是较低级与较高级对比。
• 在有序多分类Logistic回归中，假设几个二元Logistic回归中，自变量的系数相等，仅常数项不等，结果也只输出一组自变量的系数。因此，有序多分类的Logistic回归模型，必须对自变量系数相等的假设（即“比例优势”假设）进行检验（又称平行线检验）。如果不满足该假设，则考虑使用无序多分类Logistic回归。
• R包实现： MASS

dat <- read.dta("ologit.dta") #ologit.dta 为数据集
head(dat)
## one at a time, table apply, pared, and public
lapply(dat[, c("apply", "pared", "public")], table)
## three way cross tabs (xtabs) and flatten the table
ftable(xtabs(~ public + apply + pared, data = dat))
summary(dat$gpa)
sd(dat$gpa)

ggplot(dat, aes(x = apply, y = gpa)) +
  geom_boxplot(size = .75) +
  geom_jitter(alpha = .5) +
  facet_grid(pared ~ public, margins = TRUE) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

• 模型构建：

## fit ordered logit model and store results 'm'
m <- polr(apply ~ pared + public + gpa, data = dat, Hess=TRUE)
## view a summary of the model
summary(m)

• 模型评估：

## store table
(ctable <- coef(summary(m)))
## calculate and store p values
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
## combined table
(ctable <- cbind(ctable, "p value" = p))

(ci <- confint(m)) # default method gives profiled CIs
confint.default(m) # CIs assuming normality
## odds ratios
exp(coef(m))
## OR and CI
exp(cbind(OR = coef(m), ci))


sf <- function(y) {
  c('Y>=1' = qlogis(mean(y >= 1)),
    'Y>=2' = qlogis(mean(y >= 2)),
    'Y>=3' = qlogis(mean(y >= 3)))
}
(s <- with(dat, summary(as.numeric(apply) ~ pared + public + gpa, fun=sf)))

glm(I(as.numeric(apply) >= 2) ~ pared, family="binomial", data = dat)

glm(I(as.numeric(apply) >= 3) ~ pared, family="binomial", data = dat)

s[, 4] <- s[, 4] - s[, 3]
s[, 3] <- s[, 3] - s[, 3]
s # print

plot(s, which=1:3, pch=1:3, xlab='logit', main=' ', xlim=range(s[,3:4]))

newdat <- data.frame(
  pared = rep(0:1, 200),
  public = rep(0:1, each = 200),
  gpa = rep(seq(from = 1.9, to = 4, length.out = 100), 4))
newdat <- cbind(newdat, predict(m, newdat, type = "probs"))

##show first few rows
head(newdat)

lnewdat <- melt(newdat, id.vars = c("pared", "public", "gpa"),
                variable.name = "Level", value.name="Probability")
## view first few rows
head(lnewdat)

ggplot(lnewdat, aes(x = gpa, y = Probability, colour = Level)) +
  geom_line() + facet_grid(pared ~ public, labeller="label_both")