逻辑回归属于广义线性模型(Generalized linear model, GLM)
特征：
1.输入特征进行线性组合
2.输出都是指数分布一类的相关概率分布：正态分布、泊松分布、二项分布
3.输出分布和输入特征的线性组合通过函数关联
一文搞懂极大似然估计
贝叶斯公式的直观理解(先验概率/后验概率)
一篇文章完全弄懂Logistic回归（含极大似然估计详细推导和实现代码）
{高中生能看懂的}梯度下降是个啥？
逻辑回归笔记总结

1.利用线性回归进行分类

用简单线性回归进行分类，增加了两个极端点会导致结果变化。

set.seed(8967563)
##构建随机的分类数据
no_funding = rnorm(200, 15,15)
#把小于0的数字变成0，大于0的数字不变
#sapply()函数的作用是将列表、向量或数据帧作为输入，并以向量或矩阵的形式给出输出。sapply(X, FUN)：X向量或对象，FUN作用于x中每个元素的函数。
no_funding = sapply(no_funding,function(x) max(x,0)) 
#两列组合，一列是0，一列是no_funding
no_funding_df = cbind(0,no_funding)
funding = rnorm(120,50,15)
funding = sapply(funding,function(x) max(x,0))
#两列组合，一列是1，一列是funding
funding_df = cbind(1,funding)
#以上构建了两个类别，均值大的是1类，均值小的是0类
#行绑定，前面120行是funding_df，属于1类；后200行是no_funding_df，属于0类
funding_ds = rbind(funding_df, no_funding_df)
fdata = data.frame(funding_ds)
#列命名
names(fdata)=c("y","x1")
#sample随机抽取行，即打乱行的顺序
fdata = fdata[sample(nrow(fdata)),]
library(data.table)
#设置table属性，改行名
(setattr(fdata, "row.names", 1:nrow(fdata)))
head(fdata)

myfit = lm(y~x1, data=fdata)
int0.5 <- (0.5 - coef(myfit)[1]) / coef(myfit)[2] 
#coef(myfit)[1]截距，coef(myfit)[2]斜率
#用大于0.5的y是1，小于0.5的y是0来判断，找到逻辑回归的临界点 34.05678

library(ggplot2)
p <- qplot(x1,y,data=fdata) 
p <- p + ggtitle("Classifying with Simple Linear Regression")  
p <- p + geom_abline(intercept=coef(myfit)[1],slope=coef(myfit)[2],  size=1.2, show_guide=T) 
p <- p + geom_vline(xintercept = int0.5, linetype = "longdash") #中间的虚线
p <- p + theme(plot.title = element_text(lineheight=.8, face="bold"))
p

#增加两行数据，明显偏离(1,50)的均值
skewed_fdata=rbind(fdata,c(1,120),c(1,150))
myfit2 = lm(y~x1, data=skewed_fdata)
new_int0.5 <- (0.5 - coef(myfit2)[1]) / coef(myfit2)[2]
#因为两个很大的偏离点c(1,120),c(1,150)，逻辑回归的临界点变大了

p2 <- qplot(x1,y,data=skewed_fdata) 
p2 <- p2 + ggtitle("Classifying with Simple Linear Regression")  
p2 <- p2 + geom_abline(intercept=coef(myfit2)[1],slope=coef(myfit2)[2],  size=1.2, show_guide=T, linetype = "dashed") 
p2 <- p2 + geom_abline(intercept=coef(myfit)[1],slope=coef(myfit)[2],  size=1.2, show_guide=T) 
p2 <- p2 + theme(plot.title = element_text(lineheight=.8, face="bold"))
p2

直接使用线性回归模型来进行判断，鲁棒性非常差。

逻辑函数

例如，用头发长短预测性别（male/female）
假设需要知道𝑝 (𝑋) = Pr (𝑌 = 1|𝑋) 和 𝑋的关系𝑝 (𝑋) = 𝛽0 + 𝛽1𝑋
预测某些𝑋的𝑝(𝑋) > 0，或者𝑝 (𝑋)< 0

logit_seq=seq(-10,10,0.01) #-10到10，以0.01为步长，2001个数据
logit_f = 1/(1+exp(-logit_seq))

p <- qplot(logit_seq,logit_f) #设置x轴，y轴
p <- p + ggtitle("The Logistic Function")  #坐标名称
p <- p + theme(plot.title = element_text(lineheight=.8, face="bold"))
p <- p + xlab("x")  
p <- p + ylab("f(x)")  
p

y的取值范围在[0，1]之间，在远离0的地方函数的值会很快接近0或者1。当x>0时，y取值为1的可能性变大，当x<0时，y取值的0的可能性变大

2.预测心脏病

#读取数据
heart <- read.table("heart.dat", quote="\"")
names(heart) <- c("AGE", "SEX", "CHESTPAIN", "RESTBP", "CHOL", "SUGAR", "ECG", "MAXHR", "ANGINA", "DEP", "EXERCISE", "FLUOR", "THAL", "OUTPUT")

#将一些特征因子化
heart$CHESTPAIN = factor(heart$CHESTPAIN)
heart$ECG = factor(heart$ECG)
heart$THAL = factor(heart$THAL)
heart$EXERCISE = factor(heart$EXERCISE)
#将输出范围限制在0和1之间（之前是1和2）
heart$OUTPUT = heart$OUTPUT-1

#划分训练集和测试集
library(caret)
set.seed(987954)
#按照OUTPUT进行索引，提取85%作为训练集
heart_sampling_vector <- createDataPartition(heart$OUTPUT, p = 0.85, list = FALSE)
# 获取训练集数据
heart_train <- heart[heart_sampling_vector,]
heart_train_labels <- heart$OUTPUT[heart_sampling_vector]
# 获取测试集数据
heart_test <- heart[-heart_sampling_vector,]
heart_test_labels <- heart$OUTPUT[-heart_sampling_vector]

#进行逻辑回归
heart_model = glm(OUTPUT~.,data=heart_train,family=binomial("logit"))

#查看回归后的结果
summary(heart_model)

从上面的输出，可以看到：
glm函数会自动将一些因子化的变量根据level的数量自动转换为n-1个虚拟变量。其中n=level的数量。（caret包中的dummyVars函数也有同样的功能）。

注意到，年龄的系数是负的，这跟常识不符（年龄越大，得心脏病的几率越大），因此，怀疑存在多重共线性的问题。

heart_model2 = glm(OUTPUT~AGE,data=heart_train,family=binomial("logit"))
summary(heart_model2)

单独对年龄进行回归，发现年龄跟实际的心脏病是成正相关的，且P-value 变小了，显著性明显增强，因此更能证明heart_model2存在多重共线性问题。

3.评估逻辑回归模型

◾ 偏差残差
由于逻辑回归中，

◾ 空偏差
是指没有任何特征，只有常数项𝛽0的偏差。类似于线性回归中的总平方和（TSS）。其对应关系如下：

##【以下全部都是3大评价函数R方、P值、方差的准备函数】
#【准备函数1】
# 计算个体观测数据的对数似然
log_likelihoods <- function(y_labels,y_probs) {
  y_a = as.numeric(y_labels)
  y_p = as.numeric(y_probs)
  y_a * log(y_p) + (1 - y_a) * log(1 - y_p)
}

#【准备函数2】
# 计算整个数据集的对数似然
dataset_log_likelihood <- function(y_labels,y_probs) {
  sum(log_likelihoods(y_labels,y_probs))
}

#【准备函数3】
# 计算一个数据的偏差
deviances <- function(y_labels,y_probs) {
  -2*log_likelihoods(y_labels,y_probs)
}

#【准备函数4】
# 计算整个数据集的偏差，所有观测数据的偏差和
dataset_deviance <- function(y_labels, y_probs) {
  sum(deviances(y_labels,y_probs))
}

#【准备函数5】
# 计算模型的偏差
model_deviance <- function(model, data, output_column) {
  y_labels = data[[output_column]]
  y_probs = predict(model, newdata=data, type = "response")
  dataset_deviance(y_labels, y_probs)
}

#【准备函数6】
# 计算空模型偏差
null_deviance <- function(data, output_column) {
  y_labels = data[[output_column]]
  y_probs = mean(data[[output_column]])
  dataset_deviance(y_labels, y_probs)
}

#【评价方式1】Rseudo R^2,代表解释了多少比例的空偏差
model_pseudo_r_squared <- function(model, data, output_column) {
  1 - (model_deviance(model, data, output_column)/null_deviance(data, output_column))
}

model_pseudo_r_squared(heart_model,data=heart_train,output_column="OUTPUT")

#【评价方式2】P值
model_chi_squared_p_value <-  function(model, data, output_column) {
  null_df = nrow(data) - 1
  model_df = nrow(data) - length(model$coefficients)
  difference_df = null_df - model_df
  null_deviance = null_deviance(data, output_column)
  m_deviance = model_deviance(model, data, output_column)
  difference_deviance = null_deviance - m_deviance
  pchisq(difference_deviance, difference_df,lower.tail=F)
}

#【评价方式3】残差
model_deviance_residuals <- function(model, data, output_column) {
  y_labels = data[[output_column]]
  y_probs = predict(model, newdata=data, type = "response")
  residual_sign = sign(y_labels - y_probs)
  residuals = sqrt(deviances(y_labels,y_probs))
  residual_sign*residuals
}

4.模型性能测试

有了模型后，就可以来进行预测：

#对训练集进行预测，查看模型在训练集上的拟合程度：
train_predictions = predict(heart_model, newdata=heart_train, type="response")  #type=response则返回为概率值
#对于每一个观测值，大于0.5的判断为有病（1），小于等于0.5的判断为没有病
train_class_predictions = as.numeric(train_predictions>0.5)
#统计预测正确率
mean(train_class_predictions==heart_train$OUTPUT)
#查看预测错误的数据
train_predictions[train_class_predictions!=heart_train$OUTPUT]

#对测试集进行预测，查看我们模型在训练集上的拟合程度：
test_predictions = predict(heart_model, newdata=heart_test, type="response")
test_class_predictions = as.numeric(test_predictions>0.5)
mean(test_class_predictions==heart_test$OUTPUT) 

返回数值越高，意味着其精度越高。

5.分类指标

• 查全率（Recall）：正确1/真实1
• 查准率（Precise）：真实1/(真实1+错误1)
• 特异性（Specificity）：真实0/(真实0+错误0)

在这里引入一个评判分类模型好坏的统计量——AUC。
ROC（Receiver Operating Characteristic Curve），即受试者工作特征曲线。
该曲线的横坐标为假阳性率（False Positive Rate, FPR），FPR=FP/(FP+FN)，纵坐标为真阳性率（True Positive Rate, TPR）， TPR=TP/(TP+TN)，可以看到，其实真阳性率就是之前说到的查全率。
ROC的图形一般如下：

对于同一个模型，判断的阈值不同，对应的TPR和FPR也会不同。所有阈值对应的PR和FPR形成的点连起来，就形成了ROC曲线。

#计算混淆矩阵
confusion_matrix <- table(predicted = train_class_predictions, actual = heart_train$OUTPUT)
#计算查准率
precision <- confusion_matrix[2,2]/sum(confusion_matrix[2,])
#计算查全率
recall <- confusion_matrix[2,2]/sum(confusion_matrix[,2])
#计算F分数
f = 2 * precision * recall / (precision + recall)

#画ROC曲线或者P-R曲线
library(ROCR)
train_predictions = predict(heart_model, newdata=heart_train, type="response")
pred <- prediction(train_predictions, heart_train$OUTPUT)
perf <- performance(pred, measure = "prec", x.measure = "rec") #画p-R曲线
perf <- performance(pred,"tpr","fpr"); #画ROC曲线

plot(perf, main = "Precision-Recall Curve for Heart Model", lwd = 2)

一般会希望真阳性率越高越好，而假阳性率越低越好，也就是点越靠近左上角，结果是越好的。
所以，定义了曲线下面积 AUC 。这个面积越大就说明整体更靠近左上，效果也就越好。

6.正则化与过拟合

对于逻辑回归，当特征数较多时，也会有可能存在过拟合的情况。在线性回归中的有偏估计（Lasso和岭回归），对于逻辑回归同样有用：

library(glmnet)
#使用model.matrix创建特征矩阵，确保各列都是数值型
heart_train_mat <- model.matrix(OUTPUT~., heart_train)[,-1]
#给定lamda范围,从10^8 到 10^-4,平均生成250个lamda
lambdas <- 10 ^ seq(8,-4,length=250)
#lasso回归
heart_models_lasso= glmnet(heart_train_mat,heart_train$OUTPUT,alpha=1,lambda=lambdas, family="binomial")

#使用cv.glmnet选择最优的𝜆,family="binomial"二分类
lasso.cv <- cv.glmnet(heart_train_mat,heart_train$OUTPUT,alpha=1,lambda=lambdas,family="binomial")
lambda_lasso <- lasso.cv$lambda.min
lambda_lasso

#输出新数据的预测值，type参数允许选择预测的类型并提供预测值，newx代表要预测的数据
predict(heart_models_lasso, type="coefficients", s = lambda_lasso)

lasso_train_predictions <- predict(heart_models_lasso, s = lambda_lasso, newx = heart_train_mat, type="response")
lasso_train_class_predictions = as.numeric(lasso_train_predictions>0.5)
mean(lasso_train_class_predictions==heart_train$OUTPUT)

heart_test_mat <- model.matrix(OUTPUT~., heart_test)[,-1]
lasso_test_predictions <- predict(heart_models_lasso, s = lambda_lasso, newx = heart_test_mat, type="response")
lasso_test_class_predictions = as.numeric(lasso_test_predictions>0.5)
mean(lasso_test_class_predictions==heart_test$OUTPUT)

7.多元逻辑回归模型

1)无序多元逻辑回归

对于无序多元逻辑回归，可建立K-1个二元回归模型（其中k表示类别个数）。给定一个数据点，K-1个模型都会运行，概率最大的类别将会被选为预测类别。对于输入点x，分类结果为各类别的概率分别为如下公式：

在R中，无序逻辑回归可以使用glmnet::multinorm函数来预测：

#读取数据
glass <- read.csv("glass.data", header=FALSE)
#重命名属性
names(glass) <- c("id","RI","Na", "Mg", "Al", "Si", "K", "Ca", "Ba", "Fe", "Type")
#去掉第一列（ID不是特征）
glass <- glass[,-1]

#划分训练集和测试集
set.seed(4365677)
glass_sampling_vector <- createDataPartition(glass$Type, p = 0.80, list = FALSE)
glass_train <- glass[glass_sampling_vector,]
glass_test <- glass[-glass_sampling_vector,]

#使用multinom函数来进行多元逻辑回归
library(nnet)
glass_model <- multinom(Type ~ ., data = glass_train, maxit=1000)
#查看模型结果
summary(glass_model)

Type中有6种输出变量水平（1，2，3，5，6，7），所以有5组系数（2，3，5，6，7）。
其原因是，多项逻辑回归的原理是，在结果变量的类别中，会选择一个类别作为参考类别，其他各个类别的系数结果都是以参考类别作为参考（即定义为1），进一步计算及进行诠释。
在二分类的逻辑回归中，也应用到了这一点，0/1变量中，默认的0是参考类别，1是感兴趣类别。但是由于二分类结果没有第三个类别，所以不需要额外提出所谓的参考和非参考类别。但是多项逻辑回归不同，有多个类别在结果变量，需要考虑参考类别的选择。

#查看训练集的回归结果
glass_predictions <- predict(glass_model, glass_train)
mean(glass_predictions==glass_train$Type)
table(predicted=glass_predictions,actual=glass_train$Type)

#查看测试集的回归结果
glass_test_predictions <- predict(glass_model, glass_test)
mean(glass_test_predictions==glass_test$Type)
table(predicted = glass_test_predictions,actual =glass_test$Type)

• 1和2预测得较差
• 3完全被1和2混为一谈
• 6被正确区分，7只有1个错误

关于模型结果的解读，参考R-多项逻辑回归／多类别逻辑回归

2)有序多元逻辑回归

对于有序逻辑回归，如果有K个类别，需要对一个二元逻辑回归模型的输出设置K-1个阈值。该模型用决定该直线斜率的一组βi系数，但是包含K-1个截距项。

预测葡萄酒品质

winequality.white <- read.csv("winequality-white.csv", sep=";")
wine <- winequality.white
#由于wine的品质有1，2，3，4，5，6，7等多种类型，但是实际上不需要这么多类型
#划分好、中、坏三种类型，小于5的归为0，大于6的归为2，其余归为1
wine$quality <- factor(ifelse(wine$quality < 5, 0, ifelse(wine$quality > 6, 2, 1)))

#划分测试集和训练集
set.seed(7644)
wine_sampling_vector <- createDataPartition(wine$quality, p = 0.80, list = FALSE)
wine_train <- wine[wine_sampling_vector,]
wine_test <- wine[-wine_sampling_vector,]

#使用polr函数来进行有序逻辑回归
library(MASS)
wine_model <- polr(quality ~., data = wine_train, Hess=T)
#查看回归结果
summary(wine_model)

#查看训练集的回归结果
wine_predictions <- predict(wine_model, wine_train)
mean(wine_predictions==wine_train$quality)
table(predicted=wine_predictions,actual=wine_train$quality)

#查看测试集的回归结果
wine_test_predictions <- predict(wine_model, wine_test)
mean(wine_test_predictions==wine_test$quality)
table(predicted = wine_test_predictions,actual =wine_test$quality)

从上面的结果可以看到，有3个类别，产生了两个截距项。
关于有序逻辑回归的更多解释，可以参考《有序logistic回归》。

对于有序多元逻辑回归，我们也可以采用多分类器的方法来进行回归预测。

wine_model2 <- multinom(quality ~ ., data = wine_train, maxit=1000)
wine_predictions2 <- predict(wine_model2, wine_test)
mean(wine_predictions2==wine_test$quality)
table(predicted=wine_predictions2,actual=wine_test$quality)

#比较两种分类器的优劣
AIC(wine_model)
AIC(wine_model2)

#有序逻辑回归，也可以使用step函数来进行参数的选择（默认是向后剔除）
wine_model3 <- step(wine_model)
wine_test_predictions3 <- predict(wine_model3, wine_test)
mean(wine_test_predictions3==wine_test$quality)
table(predicted = wine_test_predictions3,actual =wine_test$quality)