给女朋友写的生统资料_Part17

2019-06-16 本文已影响0人城管大队哈队长

逻辑斯蒂回归

我们之前提到的线性回归是利用X来预测Y，Y是连续型的数值变量。但有时候Y并不是连续型的变量，而是一种离散型的变量。或者说，更准确来说，是一种定类数据。举个《统计学习导论》书上的例子，假设现在要通过一个急救室病人的症状来预测其患病情况。我们有三种可能的诊断：中风，服药过量，癫痫发作。我们分别用数字来表示这种诊断

中风：1
服药过量：2
癫痫发作：3

这里的诊断结果就是Y，症状就是X。我们也可以用前面线性回归来做，但这样就是默认其实是一个有序的输出。但实际上，中风，服药过量，癫痫发作这三类数据虽然我们用数字来代表，但其实并不是有顺序之分的，所以用线性回归并不适合。所以，我们就可以考虑用logistic回归来解决这种分类问题。

我之所以说Y是一种离散型的变量，更准确来说是一种定类数据。是因为还有离散型变量里面还有一类数据，叫做定序数据。其虽然是离散的，但数据之间是有程度之分的，比如轻微、严重、特别严重等等。这时候，我们既不能用线性回归，也不能用logistic回归，我们应该用的是ordinal regression（定序回归）。当然，这部分绝对不会考。。。感兴趣的可以先去看看我之前推荐过的说人话的统计学

因为通常我们遇到的一般都是二元的响应变量，即预测结果只有两类（下雨不下雨，患病不患病），包括课上讲的也是二元的，所以后面我提到的都是二元的。

在做logistic回归的时候，我们也会将我们最后的二元结果用数字来表示，一个代表1，一个代表0。我们最后预测能得到的是y=1 的概率！

logistic模型

我们来看下logistic的模型
$ln(\frac{p}{1-p}) = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdot\cdot\cdot + \beta_p X_p + \epsilon$
其中 p 代表 y = 1 的概率，x 代表了不同的自变量， $\epsilon$ 表示了误差项。与线性回归模型对比，等式右边完全相同，实际上逻辑回归模型也是广义上的线性模型。而等式的左边形式更复杂了，引入了一些非线性的变换。大家可以看到，我们这里等式左边变成了一个对数，我们常称为对数发生比（log-odd）或分对数（logit）。对数里面的 $\frac{p}{1-p}$ 就是发生比（odd），取值范围可以从0到正无穷。

然后我们估计回归参数的话，就变成
$ln(\frac{\hat{p}}{1-\hat{p}}) = \hat{\beta}_0 + \hat{\beta}_1 x_1 + \hat{\beta}_2 x_2 + \cdot\cdot\cdot + \hat{\beta}_p x_p$

有时候，也会将logistic的模型写成（我不太喜欢这么写，但有时候看助教答案是这么写的）
$\hat{y}=\hat{p}=\frac{e^{\hat{\beta}_0 + \hat{\beta}_1 x_1 + \hat{\beta}_2 x_2 + \cdot\cdot\cdot + \hat{\beta}_p x_p}}{1+e^{\hat{\beta}_0 + \hat{\beta}_1 x_1 + \hat{\beta}_2 x_2 + \cdot\cdot\cdot + \hat{\beta}_p x_p}}$

$=\frac{1}{1+e^{-({\hat{\beta}_0 + \hat{\beta}_1 x_1 + \hat{\beta}_2 x_2 + \cdot\cdot\cdot + \hat{\beta}_p x_p})}}$

这里的y就直接表示 y=1 的概率了。我感觉不太直观……

下面是logistic回归的一些tip

发生比（odd）这个概念估计是会考到的。
因为我们通常来说，会把成功写作1，失败写作0。那么发生比就是成功与失败的比例，有时候就写作 odds of success。
谁是0，谁是1的设置其实是很重要的，如果设置反了，就变成预测对立的概率了。这点我在后面例子的时候再详细地讲。
logistic回归的模型显著性和参数显著性我就不讲了，一来我不太懂╮（╯＿╰）╭，二来老师上课也没怎么提。
对于我们前面提到的逻辑斯蒂回归模型而言，系数的意思（比如 $\beta_1$ ）代表的是可以解释为在所有其他预测变量保持不变的情况下， $X_1$ 增加一个单位，对数发生比的变化为 $\beta_1$ 。

R语言的实现

用两个作业里面的题目举个例子：

题目1

数据文件“Drivers.csv”为对45名司机的调查结果，其中四个变量的含义为：

x1：表示视力状况，它是一个分类变量，1表示好，0表示有问题；
x2：年龄，数值型；
x3：驾车教育，它也是一个分类变量，1表示参加过驾车教育，0表示没有；
y：一个分类型输出变量，表示去年是否出过事故，1表示出过事故，0表示没有；

把这部分数据称为test5，已经放入Part0

# 读取并整理数据
test5 <- read.table("rawdata/homework-6.5-Drivers.csv",header = T,sep = ",")
colnames(test5) <- c("eyesight","age","eduction","accident")

> head(test5)
  eyesight age eduction accident
1        1  17        1        1
2        1  44        0        0
3        1  48        1        0
4        1  55        0        0
5        1  75        1        1
6        0  35        0        1

1.请在R语言中调用logistic回归函数，计算视力状况、年龄、驾车教育与是否发生事故的logistic回归模型，并以“odds=……”的形式写出回归公式。

R语言里面调用logistic回归函数是用glm函数。glm函数其实是拟合广义线性模型的函数，logistic回归只是其中一种。所以你要加上family=binomial，代表逻辑斯蒂回归

> test5_logistic <- glm(accident ~ .,data = test5,family = binomial())
> summary(test5_logistic)

Call:
glm(formula = accident ~ ., family = binomial(), data = test5)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5636  -0.9131  -0.7892   0.9637   1.6000  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)  
(Intercept)  0.597610   0.894831   0.668   0.5042  
eyesight    -1.496084   0.704861  -2.123   0.0338 *
age         -0.001595   0.016758  -0.095   0.9242  
eduction     0.315865   0.701093   0.451   0.6523  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 62.183  on 44  degrees of freedom
Residual deviance: 57.026  on 41  degrees of freedom
AIC: 65.026

Number of Fisher Scoring iterations: 4

然后我们的logistic回归模型就是
$ln(\frac{P_{事故}}{1-P_{事故}})=0.5976-1.496X_1-0.0016X_2+0.3159X_3$
但题目要求的是odds的形式，那么我们再改写下
$\frac{P_{事故}}{1-P_{事故}}=e^{0.5976-1.496X_1-0.0016X_2+0.3159X_3}$

2.指出（1）得到的模型中哪些因素对是否发生事故有显著性影响。如果存在对是否发生事故没有显著性影响的因素，请去除这些因素后重新计算logistic回归模型，并以“p=……”的形式写出回归公式

我们根据前面结果的p-value，就发现视力是由显著性影响的，其他因素没有显著性影响。

然后我们去掉这些因素，只留下视力这个因素，再次拟合一个logistic回归模型

> test5_logistic_reduced <- glm(accident ~ eyesight,data = test5, family = binomial())
> summary(test5_logistic_reduced)

Call:
glm(formula = accident ~ eyesight, family = binomial(), data = test5)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.4490  -0.8782  -0.8782   0.9282   1.5096  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)   0.6190     0.4688   1.320   0.1867  
eyesight     -1.3728     0.6353  -2.161   0.0307 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 62.183  on 44  degrees of freedom
Residual deviance: 57.241  on 43  degrees of freedom
AIC: 61.241

Number of Fisher Scoring iterations: 4

# 可以输出下截距，虽然前面的结果也已经输出截距了
> coefficients(test5_logistic_reduced)
(Intercept)    eyesight 
  0.6190392  -1.3728110

这回是以p的形式
$P_{事故}=\frac{e^{0.6190-1.3728X_1}}{e^{0.6190-1.3728X_1}+1}$

我们还可以比较下，这两个方程在统计学上是不是有差异的，就是跟线性回归一样用anova函数

> anova(test5_logistic,test5_logistic_reduced,test = "Chisq")
Analysis of Deviance Table

Model 1: accident ~ eyesight + age + eduction
Model 2: accident ~ eyesight
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1        41     57.026                     
2        43     57.241 -2 -0.21572   0.8978

发现两个模型是没有差异的。

3.A是一名参加过驾车教育，但视力有问题的50岁老司机；B是一名没有参加过驾车教育，但视力良好的20岁新手。现在A、B都想在某保险公司投保，但按公司规定，被保险人必须满足“明年出事故的概率不高于40%”的条件才能予以承保。请预测A、B两者明年出事故的概率，并告诉保险公司谁可以投保。

这里就是用我们构建的模型去预测结果。但有一个问题就是拿哪个模型去预测呢。答案里面是用了精简以后的（去掉了不显著的变量）模型。

我们可以都看下，反正就一行代码的事情。考试的时候，我觉得用精简的模型吧。虽然我也不太清楚啥时候去掉变量。

# 构建测试数据框
# 注意列名都必须和原来的数据框（或者说你用来构建模型的变量）列名是一致的
> data <- data.frame(eyesight=c(0,1),age=c(50,20),eduction=c(1,0))
> rownames(data) <- c("A","B")
> data
  eyesight age eduction
A        0  50        1
B        1  20        0

# 利用predict进行预测
# 添加到原来数据框会更加直观一点。
# 注意要写type="response"
> data$prob <- predict(test5_logistic,data,type="response")
> data$prob1 <- predict(test5_logistic_reduced,data,type="response")
> data
  eyesight age eduction      prob prob1
A        0  50        1 0.6971400  0.65
B        1  20        0 0.2828481  0.32

题目2

这题的数据比较多，且比较麻烦。我就不放数据了，写出来只是为了讲下参考变量和训练测试集等几个概念。

这题的大致目的就是用10个自变量去做出诊断，肿瘤是否是良性的（M = malignant（恶性的）, B = benign（良性的））。总的来说，数据集有357个良性肿瘤，212个恶性肿瘤，即共有569个数据点。

> head(test_mean,n = 1)
  diagnosis radius_mean texture_mean perimeter_mean area_mean smoothness_mean compactness_mean concavity_mean concave.points_mean symmetry_mean fractal_dimension_mean
1         M       17.99        10.38          122.8      1001          0.1184           0.2776         0.3001              0.1471        0.2419                0.07871

1. Use all mean features to construct a logistic regression model

因为原始数据集是包括了mean，var等等。这里只要求用mean部分的数据。所以我们就先提取了mean那部分数据集，然后还是用glm。

# 读取，提取数据
test <- read.table("rawdata/homework-6.6-data.csv",header = T,sep = ",")
test_mean <- test[,2:12]

# 构建logistic回归模型，formula那边我用一个.代替了所有的变量，这样写写更方便
> test_logistic <- glm(diagnosis ~ . , data = test_mean, family = binomial())
Warning message:
glm.fit: fitted probabilities numerically 0 or 1 occurred 

> summary(test_logistic)

Call:
glm(formula = diagnosis ~ ., family = binomial(), data = test_mean)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.95590  -0.14839  -0.03943   0.00429   2.91690  

Coefficients:
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)             -7.35952   12.85259  -0.573   0.5669    
radius_mean             -2.04930    3.71588  -0.551   0.5813    
texture_mean             0.38473    0.06454   5.961  2.5e-09 ***
perimeter_mean          -0.07151    0.50516  -0.142   0.8874    
area_mean                0.03980    0.01674   2.377   0.0174 *  
smoothness_mean         76.43227   31.95492   2.392   0.0168 *  
compactness_mean        -1.46242   20.34249  -0.072   0.9427    
concavity_mean           8.46870    8.12003   1.043   0.2970    
concave.points_mean     66.82176   28.52910   2.342   0.0192 *  
symmetry_mean           16.27824   10.63059   1.531   0.1257    
fractal_dimension_mean -68.33703   85.55666  -0.799   0.4244    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 751.44  on 568  degrees of freedom
Residual deviance: 146.13  on 558  degrees of freedom
AIC: 168.13

Number of Fisher Scoring iterations: 9

大家可以在输出结果里面看到这10个变量。

2. Then try to reduce the number of features from your last model, construct another regression model,and you will need to write down the equation of your logistic regression model(Tips: Logit P = α+β1X1+β2X2+..+βpXp)

我们可以把显著性的变量挑出来，再次构建一个新的logistic回归模型

> test_logistic_reduced <- glm(diagnosis ~ texture_mean + area_mean + smoothness_mean + concave.points_mean, data = test_mean, family = binomial())
Warning message:
glm.fit: fitted probabilities numerically 0 or 1 occurred 
> summary(test_logistic_reduced)

Call:
glm(formula = diagnosis ~ texture_mean + area_mean + smoothness_mean + 
    concave.points_mean, family = binomial(), data = test_mean)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.31798  -0.15623  -0.04212   0.01662   2.84201  

Coefficients:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)         -23.677816   3.882774  -6.098 1.07e-09 ***
texture_mean          0.362687   0.060544   5.990 2.09e-09 ***
area_mean             0.010342   0.002002   5.165 2.40e-07 ***
smoothness_mean      59.471304  25.965153   2.290    0.022 *  
concave.points_mean  76.571210  16.427864   4.661 3.15e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 751.44  on 568  degrees of freedom
Residual deviance: 156.44  on 564  degrees of freedom
AIC: 166.44

Number of Fisher Scoring iterations: 8

回归模型为
$ln(\frac{P_{M}}{1-P_{M}})=-23.677816+0.362687*X_1+0.010342*X_2+59.471304*X_3+76.571210*X_4$

3. Use proper test to test the difference between two models

用anova就可以了

> anova(test_logistic,test_logistic_reduced,test = "Chisq")
Analysis of Deviance Table

Model 1: diagnosis ~ radius_mean + texture_mean + perimeter_mean + area_mean + 
    smoothness_mean + compactness_mean + concavity_mean + concave.points_mean + 
    symmetry_mean + fractal_dimension_mean
Model 2: diagnosis ~ texture_mean + area_mean + smoothness_mean + concave.points_mean
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       558     146.13                     
2       564     156.44 -6   -10.31   0.1122

4. You may split the data properly, use part of them to train your regression model and use another part to make predictions. Lastly, you may try to calculate the accuracy of your model.(Tips: To split the data, you can use the first 398 rows as training data, use the last 171 rows as prediction data.The predict function return a value between 0 and 1, 0.~0.5 belong to the first class, and 0.5~1 belong to second class in binary classification problems)

# 构建训练集和测试集数据
data_train <- test[1:398,]
data_test <- test[399:569,]

# 利用训练集构建logistic回归模型
> test_train_logistic <- glm(diagnosis ~ texture_mean + area_mean + smoothness_mean + concave.points_mean, data = data_train, family = binomial())
> summary(test_train_logistic)

Call:
glm(formula = diagnosis ~ texture_mean + area_mean + smoothness_mean + 
    concave.points_mean, family = binomial(), data = data_train)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.39278  -0.14454  -0.02447   0.03635   2.60665  

Coefficients:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)         -27.47397    4.74798  -5.786 7.19e-09 ***
texture_mean          0.46244    0.08434   5.483 4.19e-08 ***
area_mean             0.01082    0.00235   4.606 4.11e-06 ***
smoothness_mean      90.11221   30.96961   2.910 0.003618 ** 
concave.points_mean  59.01212   17.51779   3.369 0.000755 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 544.93  on 397  degrees of freedom
Residual deviance: 108.30  on 393  degrees of freedom
AIC: 118.3

Number of Fisher Scoring iterations: 8

# 先把测试集里面真实的诊断结果提出来
diagnosis_pred <- data_test[,"diagnosis",drop=F]

# 用模型去预测测试集的结果，并把预测的概率结果跟真实结果放一起
diagnosis_pred$pred <- predict(test_train_logistic,data_test,type="response")

# 把>0.5的定义为M，即恶性。把<0.5定义为B，即良性
diagnosis_pred$pred_result <- ifelse(data_test$pred > 0.5, "M","B")

# 然后就拿预测的结果跟真实的结果进行比较，并计算有多少是TRUE的。就是所谓的accuracy
> mean(diagnosis_pred$diagnosis == diagnosis_pred$pred_result)
[1] 0.9064327

这样一套下来，大家可能会感觉有些奇怪。在题目1的时候，1代表出事故，0代表没出事故。然后 $ln(\frac{p}{1-p})$ 里面的p代表的是y=1的概率，即出事故的概率。这一切都很顺理成章。但是在题目2的时候，M代表恶性的，B代表良性的。那为啥 $ln(\frac{p}{1-p})$ 里面的p代表的就是M呢。或者说为啥M代表的就是1，B代表的就是0呢。

以下纯属个人观点

事实上，对于二元变量，glm会确定你的响应变量里面谁是reference level。或者说，会确定谁是那个0。那么，glm是怎么确定的呢。

如果你的响应变量本身就是0和1了，那么自然0就是reference level了。

这时候，只要响应变量是0和1，那么其就不需要是因子型的变量。数值型的0和1也是可以的

如果你的响应变量是字符串型的，那么首先，你的响应变量必须是因子型的变量才可以让glm决定谁是reference level，即谁是那个0。我们有两种方法来知道谁是reference level。

利用str或者直接输出响应变量（本来应该是test_mean做示例的，我懒得再改了……）

# 可以看到 levels 后面，排在最前面的就是reference level，即B是reference level
> str(test)
'data.frame':   569 obs. of  32 variables:
 $ id                     : int  842302 842517 84300903 84348301 84358402 843786 844359 84458202 844981 84501001 ...
 $ diagnosis              : Factor w/ 2 levels "B","M": 2 2 2 2 2 2 2 2 2 2 ...
 
> test$diagnosis
……
[532] B B M B M M B B B B B B B B B B B B B B B B B B B B B B B B B M M M M M M B
Levels: B M

利用contrast

# 清除地显示了B是0，M是1
> contrasts(test$diagnosis)
  M
B 0
M 1

如果你的响应变量既不是0和1，也不是因子型的变量，就会报错

# 把B和M变成字符串型的
> test_mean$diagnosis <- as.character(test_mean$diagnosis)
> test_logistic <- glm(diagnosis ~ . , data = test_mean, family = binomial())
Error in eval(family$initialize) : y values must be 0 <= y <= 1

当然，如果你的响应变量是数值型，但还是没有变成因子型，照样还是会报错，可以自己试试。

通过上面的讲述，我们就发现了B是reference level，即是0。而M是对立的那个level，即是1。而我们logistic输出的是y=1的概率，即y=恶性肿瘤的概率。然后也刚好对应的我们前面的

# 把>0.5的定义为M，即恶性。把<0.5定义为B，即良性
diagnosis_pred$pred_result <- ifelse(data_test$pred > 0.5, "M","B")

只有y=M的概率足够大，才定义为M。这里我们设定的足够大是0.5。你也可以认为大于0.9才算出是M，这样结果就会不太一样。

关于ifelse，可以去看下《R语言实战》第二版的p435，这里不再讲述。

我们也可以把我们的响应变量直接变成0和1

> test_mean$diagnosis <- factor(test_mean$diagnosis,levels = c("B","M"), labels = c(0,1))

那么，我们该如果更改我们的reference level呢，一种方法我觉得应该可以是像上面那样，响应变量直接变成0和1，非常直观，但个人觉得不太符合逻辑……。另一种就是下面那样

test_mean$diagnosis <- relevel(test_mean$diagnosis, ref = "M")

> test_mean$diagnosis
……
[532] B B M B M M B B B B B B B B B B B B B B B B B B B B B B B B B M M M M M M B
Levels: M B

# 现在M变了0，B变成了1
> contrasts(test_mean$diagnosis)
  B
M 0
B 1

最后强调下，glm会找谁是reference level，即谁是0。但我们最终得到的概率是对立的level的概率！

参考文章

https://stats.stackexchange.com/questions/207427/confused-with-the-reference-level-in-logistic-regression-in-r

https://stackoverflow.com/questions/17772775/change-reference-group-using-glm-with-binomial-family

https://stackoverflow.com/questions/23282048/logistic-regression-defining-reference-level-in-r

给女朋友写的生统资料_Part17

逻辑斯蒂回归

logistic模型

R语言的实现

题目1

题目2

猜你喜欢

热点阅读