In [7]:
import numpy as np
import pandas as pd
import math
from scipy.stats import chisquare
from statsmodels.stats.descriptivestats import sign_test
from statsmodels.sandbox.stats.multicomp import multipletests
import scipy
import scipy as sc
from statsmodels.stats.weightstats import *
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.stats.api as sms
import statsmodels.stats.multitest as smm
from itertools import combinations
import matplotlib.pyplot as plt
%pylab inline
In [8]:
df = pd.read_csv('credit_card_default_analysis.csv')
1. Размер кредитного лимита (LIMIT_BAL). В двух группах, тех людей, кто вернул кредит (default = 0) и тех, кто его не вернул (default = 1) проверьте гипотезы: a) о равенстве медианных значений кредитного лимита с помощью подходящей интервальной оценки b) о равенстве распределений с помощью одного из подходящих непараметрических критериев проверки равенства средних. Значимы ли полученные результаты с практической точки зрения ?
In [9]:
#1. a)
df.head()
Out[9]:
In [10]:
#1. a)
# Как видно далее на графиках, распределение отличается от нормального
pylab.figure(figsize=(12,8))
pylab.subplot(2,2,1)
stats.probplot(df[df.default == 0].LIMIT_BAL, dist="norm", plot=pylab)
pylab.subplot(2,2,2)
stats.probplot(df[df.default == 1].LIMIT_BAL, dist="norm", plot=pylab)
pylab.show()
In [11]:
plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.grid()
plt.hist(df[df.default == 0].LIMIT_BAL, color = 'lightgreen',range = (0,650000))
plt.xlabel('Balance (no default)')
plt.ylabel('Number')
plt.subplot(1,2,2)
plt.grid()
plt.hist(df[df.default == 1].LIMIT_BAL, color = 'r', range = (0,650000))
plt.xlabel('Balance (with default)')
plt.show()
print 'Avg and median Balance (no default): ', df[df.default == 0].LIMIT_BAL.mean(), df[df.default == 0].LIMIT_BAL.median(),\
'Avg and median Balance (with default): ', df[df.default == 1].LIMIT_BAL.mean(), df[df.default == 1].LIMIT_BAL.median()
Полученные гистограммы показывают очевидную разницу в объеме выданного кредита в разрезе дефолта. Хвост(no default) более тяжелый, значительно больше кредитов с большими объемами, это и приведет к тому, что медианы будут отличаться значительно.
In [12]:
#Воспользуемся бустрэпом и придем к успеху
def get_bootstrap_samples(data, n_samples):
indices = np.random.randint(0, len(data), (n_samples, len(data)))
samples = data[indices]
return samples
def stat_intervals(stat, alpha):
boundaries = np.percentile(stat, [100 * alpha / 2., 100 * (1 - alpha / 2.)])
return boundaries
no_def_lim = map(np.median, get_bootstrap_samples(df[df.default == 0].LIMIT_BAL.as_matrix(), 500))
with_def_lim = map(np.mean, get_bootstrap_samples(df[df.default == 1].LIMIT_BAL.as_matrix(), 500))
print 'no_def_lim median: ', sum(no_def_lim) / float(len(no_def_lim))
print 'with_def_lim median: ', sum(with_def_lim) / float(len(with_def_lim))
delta_median_scores = map(lambda x: x[0] - x[1], zip(no_def_lim, with_def_lim))
print 'delta ', sum(delta_median_scores) / float(len(delta_median_scores))
print "95% confidence interval for the difference between medians", stat_intervals(delta_median_scores, 0.05)
Ответ 1.а: Медианные значения отличаются существенно. Доверительный интервал для разницы не содержит 0, разница статистически значима.
Воспользуемся Ранговым критерием Манна-Уитни
In [8]:
from collections import namedtuple
MannwhitneyuResult = namedtuple('MannwhitneyuResult', ('statistic', 'pvalue'))
def tiecorrect(rankvals):
arr = np.sort(rankvals)
idx = np.nonzero(np.r_[True, arr[1:] != arr[:-1], True])[0]
cnt = np.diff(idx).astype(np.float64)
size = np.float64(arr.size)
return 1.0 if size < 2 else 1.0 - (cnt**3 - cnt).sum() / (size**3 - size)
def mannwhitneyu(x, y, use_continuity=True, alternative=None):
if alternative is None:
warnings.warn("Calling `mannwhitneyu` without specifying "
"`alternative` is deprecated.", DeprecationWarning)
x = np.asarray(x)
y = np.asarray(y)
n1 = len(x)
n2 = len(y)
ranked = scipy.stats.rankdata(np.concatenate((x, y)))
rankx = ranked[0:n1] # get the x-ranks
u1 = n1*n2 + (n1*(n1+1))/2.0 - np.sum(rankx, axis=0) # calc U for x
u2 = n1*n2 - u1 # remainder is U for y
T = tiecorrect(ranked)
if T == 0:
raise ValueError('All numbers are identical in mannwhitneyu')
sd = np.sqrt(T * n1 * n2 * (n1+n2+1) / 12.0)
meanrank = n1*n2/2.0 + 0.5 * use_continuity
if alternative is None or alternative == 'two-sided':
bigu = max(u1, u2)
elif alternative == 'less':
bigu = u1
elif alternative == 'greater':
bigu = u2
else:
raise ValueError("alternative should be None, 'less', 'greater' "
"or 'two-sided'")
z = (bigu - meanrank) / sd
if alternative is None:
# This behavior, equal to half the size of the two-sided
# p-value, is deprecated.
p = scipy.stats.distributions.norm.sf(abs(z))
elif alternative == 'two-sided':
p = 2 * scipy.stats.distributions.norm.sf(abs(z))
else:
p = scipy.stats.distributions.norm.sf(z)
u = u2
# This behavior is deprecated.
if alternative is None:
u = min(u1, u2)
return MannwhitneyuResult(u, p)
print mannwhitneyu(df[df.default == 0].LIMIT_BAL.values, df[df.default == 1].LIMIT_BAL.values)
Пол (SEX): Проверьте гипотезу о том, что гендерный состав группы людей вернувших и не вернувших кредит отличается. Хорошо, если вы предоставите несколько различных решений этой задачи (с помощью доверительного интервала и подходящего статистического критерия)
In [13]:
men = df[df.SEX == 1].default
women = df[df.SEX == 2].default
fig = plt.figure()
ax = fig.add_subplot(111)
N = 2
menMeans = [men.shape[0] - men.sum(), men.sum()]
womenMeans = [women.shape[0] - women.sum(), women.sum()]
## necessary variables
ind = np.arange(N) # the x locations for the groups
width = 0.35 # the width of the bars
## the bars
rects1 = ax.bar(ind, menMeans, width,
color='lightblue',
error_kw=dict(elinewidth=2,ecolor='red'))
rects2 = ax.bar(ind+width, womenMeans, width,
color='lightgreen',
error_kw=dict(elinewidth=2,ecolor='black'))
# axes and labels
ax.set_xlim(-width,len(ind)+width)
ax.set_ylim(0,women.shape[0] + 44)
ax.set_ylabel('Number')
ax.set_title('Default number by gender')
xTickMarks = ['No default', "With Default"]
ax.set_xticks(ind+width)
xtickNames = ax.set_xticklabels(xTickMarks)
plt.setp(xtickNames, rotation=0, fontsize=10)
## add a legend
ax.legend( (rects1[0], rects2[0]), ('Men', 'Women') )
plt.show()
print men.mean(), women.mean()
Очень интересные графики, женщины чаще берут кредит и чаще его возращают (21% женщин против 24% у мужчин, кто не вернул кредит). Проверим далее значимость этих утверждений. Воспользуемся биномиальным признаком для долей.
In [14]:
#Применим интервальную оценку методом Вилсона для каждого значения из графика выше:
from statsmodels.stats.proportion import proportion_confint
print '95%% confidence interval for %% of Men, that didnt return loans: [%f, %f]' % proportion_confint(men.sum(),
men.shape[0],
method = 'wilson')
print '95%% confidence interval for %% of WoMen, that didnt return loans: [%f, %f]' % proportion_confint(women.sum(),
women.shape[0],
method = 'wilson')
Доверительные интервалы не пересекаются, очень вероятно, что различие значимо в % возращаемых кредитов в зависимости пола. Далее определимся в значимости
In [15]:
#Воспользуемся Z-критерий для разности долей (независимые выборки)
def proportions_diff_confint_ind(sample1, sample2, alpha = 0.05):
z = scipy.stats.norm.ppf(1 - alpha / 2.)
p1 = float(sum(sample1)) / len(sample1)
p2 = float(sum(sample2)) / len(sample2)
left_boundary = (p1 - p2) - z * np.sqrt(p1 * (1 - p1)/ len(sample1) + p2 * (1 - p2)/ len(sample2))
right_boundary = (p1 - p2) + z * np.sqrt(p1 * (1 - p1)/ len(sample1) + p2 * (1 - p2)/ len(sample2))
return (left_boundary, right_boundary)
def proportions_diff_z_stat_ind(sample1, sample2):
n1 = len(sample1)
n2 = len(sample2)
p1 = float(sum(sample1)) / n1
p2 = float(sum(sample2)) / n2
P = float(p1*n1 + p2*n2) / (n1 + n2)
return (p1 - p2) / np.sqrt(P * (1 - P) * (1. / n1 + 1. / n2))
def proportions_diff_z_test(z_stat, alternative = 'two-sided'):
if alternative not in ('two-sided', 'less', 'greater'):
raise ValueError("alternative not recognized\n"
"should be 'two-sided', 'less' or 'greater'")
if alternative == 'two-sided':
return 2 * (1 - scipy.stats.norm.cdf(np.abs(z_stat)))
if alternative == 'less':
return scipy.stats.norm.cdf(z_stat)
if alternative == 'greater':
return 1 - scipy.stats.norm.cdf(z_stat)
print "95%% confidence interval for a difference between proportions: [%f, %f]" %\
proportions_diff_confint_ind(men, women)
print "p-value: ", proportions_diff_z_test(proportions_diff_z_stat_ind(men, women), 'greater')
На уровне значимости 0.05 нулевая гипотеза о независимости дефолта от пола уверенно отвергается. При этом доверительный интервал для разности долей не содержит ноль, что потверждает наши рассуждения. Данные имеют практическую значимость, от есть признак пол явно имеет смысл включать в модель.
Проверьте гипотезу о том, что образование не влияет на то, вернет ли человек долг. Предложите способ наглядного представления разницы в ожидаемых и наблюдаемых значениях количества человек вернувших и не вернувших долг. Например, составьте таблицу сопряженности, где значением ячейки была бы разность между количеством человек вернувших и не вернувших долг по каждому из значений образования. Похожи ли значения в этой таблице ? Как бы вы предложили модифицировать таблицу так, чтобы привести значения ячеек к одному масштабу не потеряв в интерпретируемости ? Наличие какого образования является наилучшим индикатором того, что человек отдаст долг ? наоборт, не отдаст долг ?
In [16]:
diff_by_educ = np.array([])
for educ in xrange(0,7):
diff = df[(df.EDUCATION == educ) & (df.default == 0)].shape[0] - df[(df.EDUCATION == educ) & (df.default == 1)].shape[0]
diff_by_educ = np.append(diff_by_educ, diff)
print educ, diff, df[(df.EDUCATION == educ) ].shape[0]
#print diff_by_educ
In [17]:
fig = plt.figure()
ax = fig.add_subplot(111)
N = 7
menMeans = diff_by_educ
## necessary variables
ind = np.arange(N) # the x locations for the groups
width = 0.35 # the width of the bars
## the bars
rects1 = ax.bar(ind, menMeans, width,
color='lightblue',
error_kw=dict(elinewidth=2,ecolor='red'))
# axes and labels
ax.set_xlim(-width,len(ind)+width)
ax.set_ylabel('Number')
xTickMarks = ['doc.', 'master', 'bach', 'school grad.', 'elemt', 'other', 'no data']
ax.set_xticks(ind+width)
xtickNames = ax.set_xticklabels(xTickMarks)
plt.setp(xtickNames, rotation=35, fontsize=10)
plt.show()
Получаем, что у нас действительно неравномерное распределение по обучению, нормализуем количество, разделяя разницу на общее кол-во людей, которое получало кредит в рассматриваемой группе. Перестроим график.
In [18]:
diff_by_educ = np.array([])
for educ in xrange(0,7):
diff = 1.*(df[(df.EDUCATION == educ) & (df.default == 0)].shape[0] ) \
/df[(df.EDUCATION == educ) ].shape[0]
diff_by_educ = np.append(diff_by_educ, diff)
print educ, diff
fig = plt.figure()
ax = fig.add_subplot(111)
N = 7
menMeans = diff_by_educ
## necessary variables
ind = np.arange(N) # the x locations for the groups
width = 0.35 # the width of the bars
## the bars
rects1 = ax.bar(ind, menMeans, width,
color='lightgreen',
error_kw=dict(elinewidth=2,ecolor='red'))
# axes and labels
ax.set_xlim(-width,len(ind)+width)
ax.set_ylabel('Number')
ax.set_title('Loan return number by education')
xTickMarks = ['doc.', 'master', 'bach', 'school grad.', 'elemt', 'other', 'no data']
ax.set_xticks(ind+width)
xtickNames = ax.set_xticklabels(xTickMarks)
threshold = 1 - df.default.mean()
ax.plot([0., 6.7], [threshold, threshold], "k--", label='Mean value')
plt.setp(xtickNames, rotation=35, fontsize=10)
plt.show()
Если отбросить 0й тип обучения, в которых слишком мало данных. То elem - наибольшая вероятность возрата долга, "школьное образование" наименьшая вероятность возрата. Далее проверим значимость
In [19]:
diff_by_educ = np.array([])
diff1 = df[(df.EDUCATION == 1) & (df.default == 1)].shape[0]
diff2 = df[(df.EDUCATION == 1) & (df.default == 0)].shape[0]
diff_by_educ = (np.append(diff1, diff2)).reshape((2,1))
for educ in xrange(2,7):
diff1 = df[(df.EDUCATION == educ) & (df.default == 1)].shape[0]
diff2 = df[(df.EDUCATION == educ) & (df.default == 0)].shape[0]
k = (np.append(diff1, diff2)).reshape((2,1))
diff_by_educ = np.concatenate((diff_by_educ, k), axis=1)
print diff_by_educ
def cramers_stat(confusion_matrix):
chi2 = scipy.stats.chi2_contingency(confusion_matrix, correction= True)[0]
n = confusion_matrix.sum()
return np.sqrt(1.*chi2 / (n*(min(confusion_matrix.shape)-1)))
print "Cramer V: ", cramers_stat(diff_by_educ), "p value: ", scipy.stats.chi2_contingency(diff_by_educ, correction= True)[1]
Таким образом нулевая гипотеза об отсутствии связи образования и уровня дефолта уверенно отвергается на уровне значимости 0.05. Значение V Крамера равно 0.07
Проверьте, как связан семейный статус с индикатором дефолта: нужно предложить меру, по которой можно измерить возможную связь этих переменных и посчитать ее значение.
Значения: 0 = отказываюсь отвечать; 1 = замужем/женат; 2 = холост; 3 = нет данных
In [35]:
diff_by_mar = np.array([])
for educ in xrange(0,4):
diff = df[(df.MARRIAGE == educ) & (df.default == 0)].shape[0] - df[(df.MARRIAGE == educ) & (df.default == 1)].shape[0]
diff_by_mar = np.append(diff_by_mar, diff)
print educ, diff, df[(df.MARRIAGE == educ) ].shape[0]
Применим аналогичный подход, как и в предыдущем признаке:
In [39]:
diff_by_mar = np.array([])
for educ in xrange(0,4):
diff = 1.*(df[(df.MARRIAGE == educ) & (df.default == 0)].shape[0] ) \
/df[(df.MARRIAGE == educ) ].shape[0]
diff_by_mar = np.append(diff_by_mar, diff)
print educ, diff
fig = plt.figure()
ax = fig.add_subplot(111)
N = 4
menMeans = diff_by_mar
## necessary variables
ind = np.arange(N) # the x locations for the groups
width = 0.35 # the width of the bars
## the bars
rects1 = ax.bar(ind, menMeans, width,
color='lightgreen',
error_kw=dict(elinewidth=2,ecolor='red'))
# axes and labels
ax.set_xlim(-width,len(ind)+width)
ax.set_ylabel('Number')
ax.set_title('Loan return number by Mariage')
xTickMarks = ['Refuse', 'Married', 'Single', 'No data']
ax.set_xticks(ind+width)
xtickNames = ax.set_xticklabels(xTickMarks)
threshold = 1 - df.default.mean()
ax.plot([0., 6.7], [threshold, threshold], "k--", label='Mean value')
plt.setp(xtickNames, rotation=35, fontsize=10)
plt.show()
In [49]:
diff_by_mar = np.array([])
diff1 = df[(df.MARRIAGE == 0) & (df.default == 1)].shape[0]
diff2 = df[(df.MARRIAGE == 0) & (df.default == 0)].shape[0]
diff_by_mar = (np.append(diff1, diff2)).reshape((2,1))
for educ in xrange(1,4):
diff1 = df[(df.MARRIAGE == educ) & (df.default == 1)].shape[0]
diff2 = df[(df.MARRIAGE == educ) & (df.default == 0)].shape[0]
k = (np.append(diff1, diff2)).reshape((2,1))
diff_by_mar = np.concatenate((diff_by_mar, k), axis=1)
def cramers_stat(confusion_matrix):
chi2 = scipy.stats.chi2_contingency(confusion_matrix, correction= True)[0]
n = confusion_matrix.sum()
print n, confusion_matrix.shape
return np.sqrt(1.*chi2 / (n*(min(confusion_matrix.shape)-1)))
print "Cramer V: ", cramers_stat(diff_by_mar), "p value: ", scipy.stats.chi2_contingency(diff_by_mar, correction= True)[1]
Таким образом нулевая гипотеза об отсутствии связи образования и уровня дефолта уверенно отвергается на уровне значимости 0.05. Значение V Крамера равно 0.034. В данном случае уровень значимости значительно меньше, но при этом сильно меньше 0.05
Возраст (AGE): Относительно двух групп людей вернувших и не вернувших кредит проверьте следующие гипотезы: a) о равенстве медианных значений возрастов людей b) о равенстве распределений с помощью одного из подходящих непараметрических критериев проверки равенства средних. Значимы ли полученные результаты с практической точки зрения ?
Применим аналогичные приемы, как и в пункте 1. Я бы применил перестановочный критерий, но мой старенький комп как то плохо строит.
In [3]:
# Как видно далее на графиках, распределение отличается от нормального
pylab.figure(figsize=(12,8))
pylab.subplot(2,2,1)
stats.probplot(df[df.default == 0].AGE, dist="norm", plot=pylab)
pylab.subplot(2,2,2)
stats.probplot(df[df.default == 1].AGE, dist="norm", plot=pylab)
pylab.show()
In [4]:
ages = np.array([])
age1 = df[(df.AGE == 18) & (df.default == 1)].shape[0]
age2 = df[(df.AGE == 18) & (df.default == 0)].shape[0]
ages = (np.append(age1, age1)).reshape((2,1))
for educ in xrange(19,81):
age1 = df[(df.AGE == educ) & (df.default == 1)].shape[0]
age2 = df[(df.AGE == educ) & (df.default == 0)].shape[0]
k = (np.append(age1, age2)).reshape((2,1))
ages = np.concatenate((ages, k), axis=1)
print ages, ages.shape
In [11]:
%pylab inline
fig = plt.figure()
ax = fig.gca(projection='3d')
X = ages[0]
Y = ages[1]
xi = linspace(min(X), max(X))
yi = linspace(min(Y), max(Y))
Z = np.arange(18,81)
zi = griddata(X, Y, Z, xi, yi, interp='linear')
X, Y = np.meshgrid(X, Y)
ax.contour(xi, yi, zi)
ax.set_xlabel('No default')
ax.set_ylabel('With default')
ax.set_zlabel('Age')
plt.show()
In [13]:
%matplotlib qt
fig = plt.figure()
ax = fig.gca(projection='3d')
X = ages[0]
Y = ages[1]
xi = linspace(min(X), max(X))
yi = linspace(min(Y), max(Y))
Z = np.arange(18,81)
zi = griddata(X, Y, Z, xi, yi, interp='linear')
X, Y = np.meshgrid(X, Y)
ax.contour(xi, yi, zi)
ax.set_xlabel('No default')
ax.set_ylabel('With default')
ax.set_zlabel('Age')
plt.show()
In [5]:
print mannwhitneyu(df[df.default == 0].AGE.values, df[df.default == 1].AGE.values)
Критерий Манна-Уитни на уровне значимости 0.05 нулевую гипотезу о равенстве средних не отклоняет (Я бы применил перестановочный критерий, но комп слабоват)
In [11]:
%pylab inline
plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.grid()
plt.hist(df[df.default == 0].AGE, color = 'lightgreen',range = (15,85))
plt.xlabel('Age (no default)')
plt.ylabel('Number')
plt.subplot(1,2,2)
plt.grid()
plt.hist(df[df.default == 1].AGE, color = 'r', range = (15,85))
plt.xlabel('Age (with default)')
plt.show()
print 'Avg and median Age (no default): ', df[df.default == 0].AGE.mean(), df[df.default == 0].AGE.median(),\
'Avg and median Age (with default): ', df[df.default == 1].AGE.mean(), df[df.default == 1].AGE.median()
Изучим медианные значения, на нашинных данных они равны, построим доверительный интервал для него:
In [14]:
#Воспользуемся бустрэпом и придем к успеху
def get_bootstrap_samples(data, n_samples):
indices = np.random.randint(0, len(data), (n_samples, len(data)))
samples = data[indices]
return samples
def stat_intervals(stat, alpha):
boundaries = np.percentile(stat, [100 * alpha / 2., 100 * (1 - alpha / 2.)])
return boundaries
no_def_age = map(np.median, get_bootstrap_samples(df[df.default == 0].AGE.as_matrix(), 500))
with_def_age = map(np.mean, get_bootstrap_samples(df[df.default == 1].AGE.as_matrix(), 500))
print 'no_def_age median: ', sum(no_def_age) / float(len(no_def_age))
print 'with_def_age median: ', sum(with_def_age) / float(len(with_def_age))
delta_median_scores = map(lambda x: x[0] - x[1], zip(no_def_age, with_def_age))
print 'delta ', sum(delta_median_scores) / float(len(delta_median_scores))
print "95% confidence interval for the difference between medians", stat_intervals(delta_median_scores, 0.05)
Медианные значения отличаются существенно. Доверительный интервал для разницы не содержит 0, разница статистически значима.
In [ ]: