import pandas as pd
import numpy as np
from statsmodels.genmod.generalized_linear_model import GLM
from statsmodels.genmod.families import Binomial
from statsmodels.stats.outliers_influence import variance_inflation_factor
import matplotlib.pyplot as plt
from sklearn.metrics import roc_curve
from sklearn.metrics import roc_auc_score
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix
from sklearn.metrics import accuracy_score
from sklearn.metrics import precision_recall_fscore_support
from sklearn.metrics import accuracy_score
from sklearn.model_selection import train_test_split
from scipy.optimize import minimize
import matplotlib
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')
df=pd.read_csv('./user_action_10days.csv')
# 파일 읽기
df=df.fillna('')
# 결측치를 공백으로 제거
df['continuous']=0
# 연속으로 들었는지 확인하는 변수를 만들기 위해 0 값을 넣어 continuous열 생성
for i in df.index:
if (len(df.iloc[i,1])!=0) & (len(df.iloc[i,2])!=0)|(len(df.iloc[i,2])!=0) & (len(df.iloc[i,3])!=0)|(len(df.iloc[i,3])!=0) & (len(df.iloc[i,4])!=0)|(len(df.iloc[i,4])!=0) & (len(df.iloc[i,5])!=0)|(len(df.iloc[i,5])!=0) & (len(df.iloc[i,6])!=0)|(len(df.iloc[i,6])!=0) & (len(df.iloc[i,7])!=0)|(len(df.iloc[i,7])!=0) & (len(df.iloc[i,8])!=0)|(len(df.iloc[i,8])!=0) & (len(df.iloc[i,9])!=0)|(len(df.iloc[i,9])!=0) & (len(df.iloc[i,10])!=0):
df['continuous'][i]=1
else : df['continuous'][i]=0
# day1 ,2 혹은 day 2,3 혹은 day 3,4 ... 연속으로 들은 접속 기록이 있으면 1 없으면 0
df.drop(df[(df['day1'] =='') & (df['day2'] == '')& (df['day3'] == '')& (df['day4'] == '')& (df['day5'] == '')& (df['day6'] == '')& (df['day7'] == '')& (df['day8'] == '')& (df['day9'] == '')& (df['day10'] == '')].index, inplace=True)
# day1 부터 day10까지 아무런 값이 없는 행 삭제
df=df.replace('','NA')
# 그후 결측치가 있는 셀을 파악하기 위해 값을 NA 로 변경
list1=[]
for i in df.index:
list1.append(df.iloc[i,1]+df.iloc[i,2]+df.iloc[i,3]+df.iloc[i,4]+df.iloc[i,5]+df.iloc[i,6]+df.iloc[i,7]+df.iloc[i,8]+df.iloc[i,9]+df.iloc[i,10])
# day 1 부터 day 10까지의 행동들을 모두 더함 (NA 값 포함)
df['sum']=list1
# 다 더한 값을 SUM 열에 추가해줌
list2=[]
for i in df.index:
list2.append(df['sum'][i].count('NA'))
# 각각 결석 횟수 (수업을 하나도 듣지 않은 날)을 파악하기 위해 NA 값을 카운트함
df['결석횟수']=list2
# 결석 횟수 열을 생성해줌
df=df.replace('NA','')
# 결석 횟수를 파악했으니 NA값을 빈칸으로 변경
all_letters=df['day1'].str.cat()
action_list= list(set(list(all_letters)))
# 각 행동별 갯수를 파악하기 위해 action list 생성
for act in action_list:
col_name=act+"_cnt"
df[col_name]=df['sum'].apply(lambda x: x.count(act))
df.iloc[:, [11,12,14,15,16,17,18,19,20,21,22,23,24]]
# 최종적으로 사용할 데이터 : 구독 여부, 연속 출석 여부, 결석 횟수, 각 행동별 카운트
X = df.iloc[:, [12,14,15,16,17,18,19,20,21,22,23,24]].values
Y = df.iloc[:, 11].values
model = GLM(Y,X,family=Binomial())
results = model.fit()
# 변수별 회귀 계수 및 p-value 확인을 위해 GLM (일반화 선형모형) 함수 사용 반응변수가 이항분포이므로 family=Binomial 사용
results.summary()
# 이제 여기서 필요한 변수를 뽑기 위해 후진제거법 사용 (모든 변수를 다 넣어서 AIC 확인 후 변수를 삭제해가며 낮은 AIC가 나올때까지 반복)
Pseudo R-squ(R제곱)
독립변수가 종속변수를 얼마나 잘 설명하는지 나타냄 (x,y 사이의 관계)
Log-Likelihood
어떤 모델과 데이터가 주어졌을 때, 이 모델이 데이터를 얼마나 잘 설명하는지를 수치화 할 수 있는 척도 (모델과 데이터 사이의 관계)
Deviance: 2[L1-L0]
조화모형과 독립모형 라이클리 차이(이탈도)
L1= 대립가설(복잡모형) 라이클리후드
L0= 귀무가설(간단모형) 라이클리후드
Deviance는 카이제곱 분포를 기각역으로 한다. 이탈도(deviance ~ x^2)
H0: 모형이 맞음 변수가 연관성 없다.
AIC: -2(log-likelihood - number of parameter)
Likelihood(우도; 가능도)란? 어떤 모델(파라미터;확률분포)과 데이터가 주어졌을 때, 이 모델이 데이터를 얼마나 잘 설명하는지? 수치화 할 수 있는 척도이다.
AIC는 로그가능도함수 값이 높으면 가산점을 주고 모형에 사용된 모수가 많으면 패널티를 주는 형식으로 만들어진 척도로 작으면 작을수록 바람직한 모형이라고 판단한다. AIC가 작을 수록 좋은 모형이라고 할 수 있다.
ex) 2개의 모수를 이용한 모형의 로그가능도함수 값이-16
4개의 모수만을 이용한 모형의 로그가능도함수 값은 -15
각 모형의 AIC는 각각 36과 38.
비록 로그가능도함수 값은 두 번째 모형이 더 컸지만 모형에 사용된 모수가 더 많아 첫 번째 모형에 비해 바람직하지 않다는 판단을 할 수 있다.
후진제거법 in R studio
#fit three models
model1 <- lm(mpg ~ disp + hp + wt + qsec, data = mtcars)
model2 <- lm(mpg ~ disp + qsec, data = mtcars)
model3 <- lm(mpg ~ disp + wt, data = mtcars)
library(AICcmodavg)
#define list of models
models <- list(model1, model2, model3)
#specify model names
mod.names <- c('disp.hp.wt.qsec', 'disp.qsec', 'disp.wt')
#calculate AIC of each model
aictab(cand.set = models, modnames = mod.names)
Model selection based on AIC:
K AIC Delta_AIC AICWt Cum.Wt LL
disp.hp.wt.qsec 6 162.43 0.00 0.83 0.83 -73.53
disp.wt 4 165.65 3.22 0.17 1.00 -78.08
disp.qsec 4 173.32 10.89 0.00 1.00 -81.92
위와 같은 방법으로 AIC값 구함.
처음 모든 변수를 선택하여 AIC값을 구한 후 위 결과 값에서 p-value가 낮은 변수부터 제거하여 AIC 값 구함.
AIC 값이 낮아지다가 다시 올라갈때까지 제거하는 후진제거법 사용
['E_cnt', 'F_cnt', 'G_cnt', 'H_cnt', 'J_cnt', 'K_cnt', 'L_cnt', 'M_cnt', 'S_cnt', 'continuous', '결석횟수']
['E_cnt', 'F_cnt', 'H_cnt', 'J_cnt', 'K_cnt', 'L_cnt', 'M_cnt', 'S_cnt', 'continuous', '결석횟수']
['E_cnt', 'F_cnt', 'J_cnt', 'K_cnt', 'L_cnt', 'M_cnt', 'S_cnt', 'continuous', '결석횟수']
['F_cnt', 'J_cnt', 'K_cnt', 'L_cnt', 'M_cnt', 'S_cnt', 'continuous', '결석횟수']
['F_cnt', 'J_cnt', 'K_cnt', 'M_cnt', 'S_cnt', 'continuous', '결석횟수'] -> aic = 3334.380
['J_cnt', 'K_cnt', 'M_cnt', 'S_cnt', 'continuous', '결석횟수']
aic가 계속 감소하다가 6번에서 다시 증가. 그러므로 5번 선택함
다중공선성
- 독립변수의 P값(P-value)이 커서 개별 인자들이 유의하지 않는 경우가 있다. 이런 경우 독립변수(X)들 간에 높은 상관관계가 있다고 의심된다.
- 독립변수들간의 상관계수를 구한다.
- 분산팽창요인(VIF, Variance Inflation Factor)을 구하여 이 값이 10을 넘는다면 보통 다중공선성의 문제가 있다.
for_vif=df.loc[:, ['F_cnt', 'J_cnt', 'K_cnt', 'M_cnt', 'S_cnt', 'continuous', '결석횟수']]
for_vif
# 위에서 후진제거법을 이용하여 찾는 변수들 중 다중공선성 확인하기
vif = pd.DataFrame()
vif["VIF Factor"] = [variance_inflation_factor(for_vif.values, i) for i in range(for_vif.shape[1])]
vif["features"] = for_vif.columns
vif
# M_cnt 와 S_cnt 사이 다중공선성이 있다고 판단됨.
for_heatmap = for_vif.corr()
matplotlib.rcParams['axes.unicode_minus'] = False
plt.rcParams['font.family'] = 'Malgun Gothic'
fig, ax = plt.subplots( figsize=(7,7) )
mask = np.zeros_like(for_heatmap, dtype=np.bool)
mask[np.triu_indices_from(mask)] = True
sns.heatmap(for_heatmap,
cmap = 'RdYlBu_r',
annot = True,
mask=mask,
linewidths=.5,
cbar_kws={"shrink": .5},
vmin = -1,vmax = 1
)
plt.show()
# 변수들간 상관관계 확인
# 다중공선성이 10이 넘는 변수 = M_cnt(19.005495) S_cnt(20.975598)
# 히트맵을 그려보고 M_cnt는 S_cnt와 강한 상관관계를 보이므로
# S_cnt만 삭제 후 M_cnt는 사용하기로 결정X = df.loc[:, ['F_cnt', 'J_cnt', 'K_cnt', 'M_cnt', 'S_cnt', 'continuous', '결석횟수']].values
Y = df.iloc[:, 11].values
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size = 0.2, random_state = 0)
classifer = LogisticRegression()
classifer.fit(X_train, Y_train)
Y_pred = classifer.predict(X_test)
accuracy_score(Y_test, Y_pred)
# 0.8451380552220888
cm=confusion_matrix(Y_test,classifer.predict(X_test))
cmtx = pd.DataFrame(
confusion_matrix(Y_test, Y_pred, labels = [1, 0]),
index = ['true:subscribe(1)', 'true:subscribe(0)'],
columns = ['pred:subscribe(1)', 'pred:subscribe(0)']
)
cmtx
precision_recall_fscore_support(Y_test,Y_pred, average='binary')
# (0.703125, 0.4972375690607735, 0.5825242718446602, None)정확도는 높으나 recall값과 f1-score가 낮다고 판단됨.
cut-off들 확인.
prob1=classifer.predict_proba(X_test)
def objective_func(cutoff):
list1=[]
for i in range(len(prob1)):
if prob1[i][1]>cutoff:
list1.append(1)
else:
list1.append(0)
return 1-(accuracy_score(Y_test, list1))
res = minimize(objective_func, 1, method='nelder-mead',
options={'xatol': 1e-8, 'disp': True})
import numpy as np
cutoff_points = np.linspace(0.4, 0.6, 1000)
accuracy_inversed_results = [objective_func(cutoff) for cutoff in cutoff_points]
plt.scatter(x=cutoff_points, y=accuracy_inversed_results)
plt.xlabel("cutoff_point")
plt.ylabel("1-accuracy")정확도와 recall을 둘 다 높이는 cut off 를 찾기 위한 리스트 출력
cut off의 0.0001 단위로 정확도 값 찾음
for i in range(10000,0,-1):
a=i*0.0001
list1=[]
list2=[]
for j in range(len(prob1)):
if prob1[j][1]>a:
list1.append(1)
else:
list1.append(0)
print(accuracy_score(Y_test, list1),a)
# 각 cut-off별 accuracy 값 파악
# cut-off는 1부터 0.0001 단위로 0까지 하락
# 왼쪽 accuracy 오른쪽 cut-off 값정확도가 0.84 보다 높은 cut off 값 출력
for i in range(100000,0,-1):
a=i*0.00001
list1=[]
list2=[]
for j in range(len(prob1)):
if prob1[j][1]>a:
list1.append(1)
else:
list1.append(0)
b=accuracy_score(Y_test, list1)
if format(b, '.4f')>'0.8400':
print(a)
else: pass
# accuracy가 0.84가 넘는 cut-off값 출력Recall값 출력
for i in range(1000,0,-1):
a=i*0.001
list1=[]
list2=[]
for j in range(len(prob1)):
if prob1[j][1]>a:
list1.append(1)
else:
list1.append(0)
cm=confusion_matrix(Y_test,list1)
Recall=cm[1,1]/(cm[1,1]+cm[1,0])
print(Recall,a)
# 왼쪽 Recall 오른쪽 cut-off 값
# Recall 값은 cut-off값과 반비례 (cut-off 값이 낮아지면 recall 높아짐)Recall값과 accuracy 값이 반비례 하므로 f1 score를 높이는 것을 목적으로 함.
for i in range(1000,0,-1):
a=i*0.001
list1=[]
list2=[]
for j in range(len(prob1)):
if prob1[j][1]>a:
list1.append(1)
else:
list1.append(0)
cm=confusion_matrix(Y_test,list1)
Recall=cm[1,1]/(cm[1,1]+cm[1,0])
Precision=cm[1,1]/(cm[1,1]+cm[0,1])
f1= 2*(Precision*Recall)/(Precision+Recall)
print(f1,a)
# 왼쪽 f1 score 값 오른쪽 cut-off 값
# 위에서 accuracy가 0.84가 넘고, cut-off값 중 f1 score가 가장 높은 값 : 0.443
최적의 cut off 값 0.443으로 설정
list1=[]
for i in range(len(prob1)):
if prob1[i][1]>0.443:
list1.append(1)
else:
list1.append(0)
accuracy_score(Y_test, list1)
#0.8427370948379351
precision_recall_fscore_support(Y_test, list1, average='binary')
# (0.6602564102564102, 0.569060773480663, 0.6112759643916914, None)정리
유료 구독 전환 여부 판단을 위해 사용한 변수
- F_cnt, J_cnt, K_cnt, M_cnt, S_cnt, continuous, 결석횟수
사용 모델
- LogisticRegression
cut off값
- 0.443
accuracy_score, precision, recall, f1score
- 0.843, 0.660, 0.569, 0.611
한걸음씩 배워나갑니다
윤섭님의 열정적인 개발 여정을 응원 드립니다 ~! 새로 안내드린 블로그챌린지 규칙 확인하셔서 태그 추가해주시면 저 제니아는 더 큰 감동 받을 것 같아요...(ﾉ◕ヮ◕)ﾉ*:･ﾟ✧ 남은 2월도 다가오는 3월도 화이팅입니다 ~!