邏輯回歸屬于廣義線性模型(Generalized linear model, GLM)
特征：
1.輸入特征進行線性組合
2.輸出都是指數(shù)分布一類的相關(guān)概率分布：正態(tài)分布、泊松分布、二項分布
3.輸出分布和輸入特征的線性組合通過函數(shù)關(guān)聯(lián)
一文搞懂極大似然估計
貝葉斯公式的直觀理解(先驗概率/后驗概率)
一篇文章完全弄懂Logistic回歸（含極大似然估計詳細推導(dǎo)和實現(xiàn)代碼）
{高中生能看懂的}梯度下降是個啥？
邏輯回歸筆記總結(jié)

1.利用線性回歸進行分類

用簡單線性回歸進行分類，增加了兩個極端點會導(dǎo)致結(jié)果變化。

set.seed(8967563)
##構(gòu)建隨機的分類數(shù)據(jù)
no_funding = rnorm(200, 15,15)
#把小于0的數(shù)字變成0，大于0的數(shù)字不變
#sapply()函數(shù)的作用是將列表、向量或數(shù)據(jù)幀作為輸入，并以向量或矩陣的形式給出輸出。sapply(X, FUN)：X向量或?qū)ο?，F(xiàn)UN作用于x中每個元素的函數(shù)。
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)
#以上構(gòu)建了兩個類別，均值大的是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)
#設(shè)置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

#增加兩行數(shù)據(jù)，明顯偏離(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

直接使用線性回歸模型來進行判斷，魯棒性非常差。

邏輯函數(shù)

例如，用頭發(fā)長短預(yù)測性別（male/female）
假設(shè)需要知道?? (??) = Pr (?? = 1|??) 和 ??的關(guān)系?? (??) = ??0 + ??1??
預(yù)測某些??的??(??) > 0，或者?? (??)< 0

logit_seq=seq(-10,10,0.01) #-10到10，以0.01為步長，2001個數(shù)據(jù)
logit_f = 1/(1+exp(-logit_seq))

p <- qplot(logit_seq,logit_f) #設(shè)置x軸，y軸
p <- p + ggtitle("The Logistic Function")  #坐標(biāo)名稱
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的地方函數(shù)的值會很快接近0或者1。當(dāng)x>0時，y取值為1的可能性變大，當(dāng)x<0時，y取值的0的可能性變大

2.預(yù)測心臟病

#讀取數(shù)據(jù)
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

#劃分訓(xùn)練集和測試集
library(caret)
set.seed(987954)
#按照OUTPUT進行索引，提取85%作為訓(xùn)練集
heart_sampling_vector <- createDataPartition(heart$OUTPUT, p = 0.85, list = FALSE)
# 獲取訓(xùn)練集數(shù)據(jù)
heart_train <- heart[heart_sampling_vector,]
heart_train_labels <- heart$OUTPUT[heart_sampling_vector]
# 獲取測試集數(shù)據(jù)
heart_test <- heart[-heart_sampling_vector,]
heart_test_labels <- heart$OUTPUT[-heart_sampling_vector]

#進行邏輯回歸
heart_model = glm(OUTPUT~.,data=heart_train,family=binomial("logit"))

#查看回歸后的結(jié)果
summary(heart_model)

從上面的輸出，可以看到：
glm函數(shù)會自動將一些因子化的變量根據(jù)level的數(shù)量自動轉(zhuǎn)換為n-1個虛擬變量。其中n=level的數(shù)量。（caret包中的dummyVars函數(shù)也有同樣的功能）。

注意到，年齡的系數(shù)是負的，這跟常識不符（年齡越大，得心臟病的幾率越大），因此，懷疑存在多重共線性的問題。

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

單獨對年齡進行回歸，發(fā)現(xiàn)年齡跟實際的心臟病是成正相關(guān)的，且P-value 變小了，顯著性明顯增強，因此更能證明heart_model2存在多重共線性問題。

3.評估邏輯回歸模型

? 偏差殘差
由于邏輯回歸中，

? 空偏差
是指沒有任何特征，只有常數(shù)項??0的偏差。類似于線性回歸中的總平方和（TSS）。其對應(yīng)關(guān)系如下：

##【以下全部都是3大評價函數(shù)R方、P值、方差的準(zhǔn)備函數(shù)】
#【準(zhǔn)備函數(shù)1】
# 計算個體觀測數(shù)據(jù)的對數(shù)似然
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)
}

#【準(zhǔn)備函數(shù)2】
# 計算整個數(shù)據(jù)集的對數(shù)似然
dataset_log_likelihood <- function(y_labels,y_probs) {
  sum(log_likelihoods(y_labels,y_probs))
}

#【準(zhǔn)備函數(shù)3】
# 計算一個數(shù)據(jù)的偏差
deviances <- function(y_labels,y_probs) {
  -2*log_likelihoods(y_labels,y_probs)
}

#【準(zhǔn)備函數(shù)4】
# 計算整個數(shù)據(jù)集的偏差，所有觀測數(shù)據(jù)的偏差和
dataset_deviance <- function(y_labels, y_probs) {
  sum(deviances(y_labels,y_probs))
}

#【準(zhǔn)備函數(shù)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)
}

#【準(zhǔn)備函數(shù)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.模型性能測試

有了模型后，就可以來進行預(yù)測：

#對訓(xùn)練集進行預(yù)測，查看模型在訓(xùn)練集上的擬合程度：
train_predictions = predict(heart_model, newdata=heart_train, type="response")  #type=response則返回為概率值
#對于每一個觀測值，大于0.5的判斷為有?。?），小于等于0.5的判斷為沒有病
train_class_predictions = as.numeric(train_predictions>0.5)
#統(tǒng)計預(yù)測正確率
mean(train_class_predictions==heart_train$OUTPUT)
#查看預(yù)測錯誤的數(shù)據(jù)
train_predictions[train_class_predictions!=heart_train$OUTPUT]

#對測試集進行預(yù)測，查看我們模型在訓(xùn)練集上的擬合程度：
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) 

返回數(shù)值越高，意味著其精度越高。

5.分類指標(biāo)

• 查全率（Recall）：正確1/真實1
• 查準(zhǔn)率（Precise）：真實1/(真實1+錯誤1)
• 特異性（Specificity）：真實0/(真實0+錯誤0)

在這里引入一個評判分類模型好壞的統(tǒng)計量——AUC。
ROC（Receiver Operating Characteristic Curve），即受試者工作特征曲線。
該曲線的橫坐標(biāo)為假陽性率（False Positive Rate, FPR），F(xiàn)PR=FP/(FP+FN)，縱坐標(biāo)為真陽性率（True Positive Rate, TPR）， TPR=TP/(TP+TN)，可以看到，其實真陽性率就是之前說到的查全率。
ROC的圖形一般如下：

對于同一個模型，判斷的閾值不同，對應(yīng)的TPR和FPR也會不同。所有閾值對應(yīng)的PR和FPR形成的點連起來，就形成了ROC曲線。

#計算混淆矩陣
confusion_matrix <- table(predicted = train_class_predictions, actual = heart_train$OUTPUT)
#計算查準(zhǔn)率
precision <- confusion_matrix[2,2]/sum(confusion_matrix[2,])
#計算查全率
recall <- confusion_matrix[2,2]/sum(confusion_matrix[,2])
#計算F分?jǐn)?shù)
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)

一般會希望真陽性率越高越好，而假陽性率越低越好，也就是點越靠近左上角，結(jié)果是越好的。
所以，定義了曲線下面積 AUC 。這個面積越大就說明整體更靠近左上，效果也就越好。

6.正則化與過擬合

對于邏輯回歸，當(dāng)特征數(shù)較多時，也會有可能存在過擬合的情況。在線性回歸中的有偏估計（Lasso和嶺回歸），對于邏輯回歸同樣有用：

library(glmnet)
#使用model.matrix創(chuàng)建特征矩陣，確保各列都是數(shù)值型
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選擇最優(yōu)的??,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

#輸出新數(shù)據(jù)的預(yù)測值，type參數(shù)允許選擇預(yù)測的類型并提供預(yù)測值，newx代表要預(yù)測的數(shù)據(jù)
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表示類別個數(shù)）。給定一個數(shù)據(jù)點，K-1個模型都會運行，概率最大的類別將會被選為預(yù)測類別。對于輸入點x，分類結(jié)果為各類別的概率分別為如下公式：

在R中，無序邏輯回歸可以使用glmnet::multinorm函數(shù)來預(yù)測：

#讀取數(shù)據(jù)
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]

#劃分訓(xùn)練集和測試集
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函數(shù)來進行多元邏輯回歸
library(nnet)
glass_model <- multinom(Type ~ ., data = glass_train, maxit=1000)
#查看模型結(jié)果
summary(glass_model)

Type中有6種輸出變量水平（1，2，3，5，6，7），所以有5組系數(shù)（2，3，5，6，7）。
其原因是，多項邏輯回歸的原理是，在結(jié)果變量的類別中，會選擇一個類別作為參考類別，其他各個類別的系數(shù)結(jié)果都是以參考類別作為參考（即定義為1），進一步計算及進行詮釋。
在二分類的邏輯回歸中，也應(yīng)用到了這一點，0/1變量中，默認的0是參考類別，1是感興趣類別。但是由于二分類結(jié)果沒有第三個類別，所以不需要額外提出所謂的參考和非參考類別。但是多項邏輯回歸不同，有多個類別在結(jié)果變量，需要考慮參考類別的選擇。

#查看訓(xùn)練集的回歸結(jié)果
glass_predictions <- predict(glass_model, glass_train)
mean(glass_predictions==glass_train$Type)
table(predicted=glass_predictions,actual=glass_train$Type)

#查看測試集的回歸結(jié)果
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預(yù)測得較差
• 3完全被1和2混為一談
• 6被正確區(qū)分，7只有1個錯誤

關(guān)于模型結(jié)果的解讀，參考R-多項邏輯回歸／多類別邏輯回歸

2)有序多元邏輯回歸

對于有序邏輯回歸，如果有K個類別，需要對一個二元邏輯回歸模型的輸出設(shè)置K-1個閾值。該模型用決定該直線斜率的一組βi系數(shù)，但是包含K-1個截距項。

預(yù)測葡萄酒品質(zhì)

winequality.white <- read.csv("winequality-white.csv", sep=";")
wine <- winequality.white
#由于wine的品質(zhì)有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)))

#劃分測試集和訓(xùn)練集
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函數(shù)來進行有序邏輯回歸
library(MASS)
wine_model <- polr(quality ~., data = wine_train, Hess=T)
#查看回歸結(jié)果
summary(wine_model)

#查看訓(xùn)練集的回歸結(jié)果
wine_predictions <- predict(wine_model, wine_train)
mean(wine_predictions==wine_train$quality)
table(predicted=wine_predictions,actual=wine_train$quality)

#查看測試集的回歸結(jié)果
wine_test_predictions <- predict(wine_model, wine_test)
mean(wine_test_predictions==wine_test$quality)
table(predicted = wine_test_predictions,actual =wine_test$quality)

從上面的結(jié)果可以看到，有3個類別，產(chǎn)生了兩個截距項。
關(guān)于有序邏輯回歸的更多解釋，可以參考《有序logistic回歸》。

對于有序多元邏輯回歸，我們也可以采用多分類器的方法來進行回歸預(yù)測。

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)

#比較兩種分類器的優(yōu)劣
AIC(wine_model)
AIC(wine_model2)

#有序邏輯回歸，也可以使用step函數(shù)來進行參數(shù)的選擇（默認是向后剔除）
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)