2020/01/26
RからPythonへの道(10)
前回の重回帰分析に続き、今回は「9. 回帰分析(ロジスティック回帰1)」について、RとPythonで計算して行きたいと思います。今回は書籍『Rによるデータマイニング入門』を教材に利用しました。 ロジスティック回帰は非線形回帰の一つで以下の式で表されます。
データはUCI Machine Learning RepositoryからSpambaseデータをダウンロードして使いました。メールの特徴からスパムメールか否かの二値判断をするものです。まずは、Rのコードです。データをプロットして比較のための単回帰分析を行い、フィッティングの悪いことを確認した上で、ロジスティック回帰を行います。
# Logistic Regression生データのみ、単回帰フィッティング、ロジスティック回帰フィッティングのグラフは以下の通りです。
library(ggplot2)
library(dplyr)
# datafile download
DLfile = tempfile()
url = "https://archive.ics.uci.edu/ml/machine-learning-databases/spambase/spambase.zip"
download.file(url, DLfile)
unzip(DLfile)
# data read
spambase = read.csv("spambase.data", header = FALSE)
colnames(spambase) = read.table("spambase.names", skip = 33, sep=":", comment.char = "")[,1]
colnames(spambase)[ncol(spambase)]="spam"
dim(spambase)
# basic statistics
spambase %>% group_by(spam) %>% summarise(count=n(), med=median(word_freq_your), mean=mean(word_freq_your), sd=sd(word_freq_your))
# graph plot(Raw data)
qplot(word_freq_your, spam, data=spambase, alpha=I(0.03))
# single regression
spam.your.lm = lm(spam~word_freq_your, data=spambase)
summary(spam.your.lm)
qplot(word_freq_your, spam, data = spambase, alpha=I(0.03), xlim = c(0, 10), ylim = c(0,1.5))+geom_smooth(method = "lm")
# logistic regression
spam.your.lst = glm(spam~word_freq_your, data=spambase, family="binomial")
summary(spam.your.lst)
a = coef(spam.your.lst)[2]
b = coef(spam.your.lst)[1]
lst.fun=function(x){
1/(1+exp(-(a*x+b)))
}
qplot(word_freq_your, spam, data = spambase, alpha=I(0.03), xlim = c(-5, 15))+stat_function(fun = lst.fun, geom="line",color="blue")
また、数値出力結果は以下の通りです。> dim(spambase)単回帰分析の調整済み決定係数が0.1467でフィッティングが良くないことが分かります。
[1] 4601 58
>
> # basic statistics
> spambase %>% group_by(spam) %>% summarise(count=n(), med=median(word_freq_your), mean=mean(word_freq_your), sd=sd(word_freq_your))
# A tibble: 2 x 5
spam count med mean sd
<int> <int> <dbl> <dbl> <dbl>
1 0 2788 0 0.439 1.03
2 1 1813 1.19 1.38 1.23
> # single regression
> summary(spam.your.lm)
Call:
lm(formula = spam ~ word_freq_your, data = spambase)
Residuals:
Min 1Q Median 3Q Max
-1.9382 -0.2677 -0.2677 0.4999 0.7322
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.267750 0.008027 33.35 <2e-16 ***
word_freq_your 0.155966 0.005543 28.14 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4514 on 4599 degrees of freedom
Multiple R-squared: 0.1469, Adjusted R-squared: 0.1467
F-statistic: 791.7 on 1 and 4599 DF, p-value: < 2.2e-16
> # logistic regression
> summary(spam.your.lst)
Call:
glm(formula = spam ~ word_freq_your, family = "binomial", data = spambase)
Deviance Residuals:
Min 1Q Median 3Q Max
-4.0211 -0.7549 -0.7549 1.1064 1.6701
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.10971 0.04235 -26.20 <2e-16 ***
word_freq_your 0.85845 0.03682 23.32 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 6170.2 on 4600 degrees of freedom
Residual deviance: 5407.3 on 4599 degrees of freedom
AIC: 5411.3
Number of Fisher Scoring iterations: 4
次に、Pythonのコードです。
# Logistic Regression生データのみ、単回帰フィッティング、ロジスティック回帰フィッティングのグラフは以下の通りです。当然ですが、Rの結果と同じグラフが得られました。
import requests, zipfile
from io import StringIO
import io
import pandas as pd
from sklearn import linear_model
import matplotlib.pyplot as plt
import numpy as np
import math
# data file download
url = 'https://archive.ics.uci.edu/ml/machine-learning-databases/spambase/spambase.zip'
r = requests.get(url, stream = True)
z = zipfile.ZipFile(io.BytesIO(r.content))
z.extractall()
# data read
spambase = pd.read_csv("spambase.data", header=None)
spambase_colnames = []
line_temp = []
with open("spambase.names","r") as f:
line_temp = f.readlines()
for i in range(33,len(line_temp)):
spambase_colnames.append(line_temp[i].split(':')[0])
spambase_colnames.append("spam")
spambase.columns = spambase_colnames
print(spambase.shape)
# basic statistics
print("count")
print(spambase.groupby('spam')['word_freq_your'].count())
print("median")
print(spambase.groupby('spam')['word_freq_your'].median())
print("mean")
print(spambase.groupby('spam')['word_freq_your'].mean())
print("std")
print(spambase.groupby('spam')['word_freq_your'].std())
# data plot
X = spambase.loc[:,['word_freq_your']].values
Y = spambase.loc[:,['spam']].values
plt.scatter(X, Y, alpha=0.03)
plt.show()
# Linear Regression (for comparison)
reg = linear_model.LinearRegression()
reg.fit(X, Y)
print('Slope: ', reg.coef_)
print('Intercept:', reg.intercept_)
print('Coefficient of determination: ',reg.score(X, Y))
plt.xlim(0.0, 10.0)
plt.ylim(0.00, 1.50)
plt.scatter(X, Y, alpha=0.03)
plt.xlabel('word_freq_your')
plt.ylabel('spam')
plt.plot(X, reg.predict(X), color= 'red') # fitting line
plt.show()
# Logistic Regression
reg1 = linear_model.LogisticRegression()
reg1.fit(X, Y)
print('Slope: ', reg1.coef_)
print('Intercept:', reg1.intercept_)
a = float(reg1.coef_)
b = float(reg1.intercept_)
def LogisticRegressionPredict(x):
return 1.0 /(1 + math.exp(-(a*x+b)))
XX = np.linspace(-5, 15,1001)
YY = np.array([LogisticRegressionPredict(i) for i in XX])
plt.xlim(-5.0, 15.0)
plt.ylim(0.00, 1.00)
plt.scatter(X, Y, alpha=0.03)
plt.xlabel('word_freq_your')
plt.ylabel('spam')
plt.plot(XX, YY, linestyle = 'None', marker=".", color= 'red') # fitting line
plt.show()
また、数値出力結果は以下の通りです。#data sizeこれも、Rの結果と同じです。ロジスティック回帰は今まで使ったことがなかったので、大変勉強になりました。RでもPythonでもやり方を理解すればすぐにでも使えそうですね。
(4601, 58)
count
spam
0 2788
1 1813
Name: word_freq_your, dtype: int64
median
spam
0 0.00
1 1.19
Name: word_freq_your, dtype: float64
mean
spam
0 0.438702
1 1.380370
Name: word_freq_your, dtype: float64
std
spam
0 1.025167
1 1.227385
Name: word_freq_your, dtype: float64
# Linear Regression (for comparison)
Slope: [[0.15596597]]
Intercept: [0.26774962]
Coefficient of determination: 0.14686816024270333
# Logistic Regression
Slope: [[0.85729313]]
Intercept: [-1.10887145]
『RからPythonへの道』バックナンバー
(1) はじめに
(2) 0. 実行環境(作業環境)
(3) 1. PythonからRを使う方法
2. RからPythonを使う方法
(4) 3. データフレーム
(5) 4. ggplot
(6) 5.行列
(7) 6.基本統計量
(8) 7. 回帰分析(単回帰)
(9) 8. 回帰分析(重回帰)
_200124.png/)
_200124.png/)
_200124.png/)
_200124.png/)
