【R语言】多项式逻辑回归logistic

require(foreign)
require(nnet)
require(ggplot2)
require(reshape2)

# 失效下载本地数据
ml <- read.dta("https://stats.idre.ucla.edu/stat/data/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")
# final negative log-likelihood 179.981726
test <- multinom(prog2 ~ ses + write, data = ml)
summary(test)
# 系数除以标准误等于t值
z <- summary(test)$coefficients/summary(test)$standard.errors;z
# 2-tailed z test
p <- (1 - pnorm(abs(z), 0, 1)) * 2;p
# extract the coefficients from the model and exponentiate
# 从模型中提取系数并求指数相当于OR值
exp(coef(test))


#  【predicted probabilities】
# fitted（拟合）是在*给定样本上*做预测，
#  而predict（预测）是在*新的样本上*做预测。

# 使用预测的概率来帮助您理解模型{1}
head(pp <- fitted(test))
# 创建一个变量变化的小数据集，同时保持另一个常数
dses <- data.frame(ses = c("low", "middle", "high"), 
                   write = mean(ml$write))
predict(test, newdata = dses, "probs")

# 使用预测的概率来帮助您理解模型{2}
# 使用预测概率理解模型的另一种方法是
# 查看每个Ses级别内连续预测变量写入的不同值的平均预测概率
summary(ml$write)
dwrite <- data.frame(ses = rep(c("low", "middle", "high"), each = 41), 
                     write = rep(c(30:70),3))
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)

# 宽数据转长
lpp <- melt(pp.write, id.vars = c("ses", "write"), 
            value.name = "probability")
head(lpp) 
# 作图看概率趋势
ggplot(lpp, aes(x = write, y = probability, colour = ses)) + 
  geom_line() + facet_grid(variable ~., scales = "free")